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I.  INRT ODCUCTION 

This  Idea  Award  (PC040282,  entitled  “Prostate  Dose  Escalation  by  Innovative  Inverse  Planning- 
Driven  IMRT”)  was  awarded  to  the  principal  investigator  (PI)  for  the  period  of  Nov  1,  2004 — Oct.  31, 
2007.  This  is  the  annual  report  for  the  first  funding  period  (Nov.  1,  2004  -  Oct.  31,  2005).  The  goal  of  this 
project  is  to  improve  current  prostate  IMRT  by  establishing  a  novel  inverse  planning  framework  with  non- 
uniform  intra-structural  penalty  distribution.  The  specific  aims  of  this  project  are  (1)  To  establish  an 
effective  inverse  planning  algorithm  with  spatially  non-uniform  importance  factors  for  prostate  IMRT;  and 
(2)  To  demonstrate  the  impact  of  the  novel  inverse  planning  formalism  for  prostate  irradiation  by  using  15 
previously  treated  prostate  cases.  Under  the  generous  support  from  the  U.S.  Army  Medical  Research  and 
Materiel  Command  (AMRMC),  the  PI  has  assembled  a  rigorous  research  team  and  setup  a  necessary 
infrastructure  needed  for  the  proposed  research.  We  have  contributed  significantly  to  prostate  cancer 
research  by  applying  physics  and  engineering  knowledge  to  prostate  cancer  research.  A  number  of 
significant  conference  abstracts  and  refereed  papers  have  been  resulted  from  the  support.  The  preliminary 
data  obtained  under  the  support  of  the  grant  has  also  enabled  the  PI  to  start  new  research  initiatives  and 
significantly  advanced  the  prostate  radiation  therapy.  In  this  report,  the  past  year’s  research  activities  of  the 
PI  are  highlighted. 


II. RESEARCH  AND  ACCOMPLISHMENTS 

The  success  of  prostate  Intensity  modulated  radiation  therapy  (IMRT)  depends  critically  on  the 
performance  of  inverse  planning  system.  As  it  is  practiced  now,  however,  the  capacity  of  IMRT  is  greatly 
underutilized  because  of  the  inferior  performance  of  the  inverse  planning  algorithms.  While  many  techniques 
exist  for  inverse  planning  and,  in  all  cases,  optimizations  are  claimed  to  be  successful,  the  plans  computed  by 
these  so-called  optimization  systems  are  often  out  of  the  expectation  of  the  physicians  and  multiple  trial-and- 
errors  are  required.  On  a  more  fundamental  level,  we  recently  found  that  the  existing  inverse  planning  suffers 
seriously  from  the  tacit  ignorance  of  intra-structural  tradeoff.  As  a  result  of  this  deficiency,  the  IMRT  plans 
generated  by  these  systems  for  prostate  treatment  are,  at  best,  sub-optimal  and  our  endeavor  of  providing  the 
best  possible  patient  care  is  compromised.  We  are  systematically  investigating  the  role  of  this  unexplored 
issue  in  inverse  planning  and  develop  innovative  techniques  to  incorporate  the  intra-structural  tradeoff  into 
prostate  IMRT  inverse  planning.  By  adequately  modeling  the  effect,  it  is  anticipated  that  the  dose  to  the 
prostate  tumor  volume  can  be  escalated  by  more  than  10%  while  maintaining  the  radiation  toxicity  at  the 
current  level  of  IMRT  prostate  treatment  (conversely,  the  dose  to  the  rectum  and/or  bladder  can  be  reduced 
greatly  if  the  target  dose  is  kept  at  the  present  level).  If  successful,  the  conformality  of  IMRT  prostate 
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treatment  will  no  longer  be  set  by  the  non-optimal  performance  of  the  dose  optimization  algorithms,  but  by 
the  physical  limit  of  IMRT.  Toward  the  general  goal  of  improving  IMRT  inverse  planning  and  treatment  of 
prostate  cancer,  we  have  done  substantial  work  in  the  past  year.  The  research  is  sorted  into  four  categories 
and  summarized  below. 

Using  a  priori  clinical  outcome  data  to  improve  inverse  planning.  IMRT  plan  ranking  is,  to  a  certain 
degree,  subjective  and  the  solution  depends  strongly  the  selection  of  models  for  inverse  planning.  Voxels 
within  a  target  or  a  sensitive  structure  volume  are  generally  not  equivalent  in  achieving  their  dosimetric 
goals  in  IMRT  planning.  Inverse  planning  objective  function  should  not  only  balance  the  competing 
objectives  of  different  structures  but  also  that  of  the  individual  voxels  in  various  structures.  While  it  is 
permissible  for  each  voxel  to  have  its  own  importance  value,  a  challenging  problem  is  how  to  obtain  a 
sensible  set  of  importance  factors  with  a  manageable  amount  of  computing.  One  of  the  most  effective  ways 
to  accomplish  a  spatially  non-uniform  penalty  scheme  is  to  incorporate  existing  clinical  knowledge  into 
inverse  planning  process.  We  have  developed  a  formalism  to  integrate  clinical  endpoint  data  to  better  guide 
the  inverse  planning  optimization  process.  In  this  work,  we  have  developed  a  new  formalism  that  is 
biologically  more  sensible  yet  clinically  practical  for  IMRT  inverse  planning.  In  this  formalism  the  dose- 
volume  status  of  a  structure  is  characterized  by  using  the  concept  of  effective  volume  in  the  voxel  domain 
and  an  objective  function  with  incorporation  of  the  volumetric  information  is  then  constructed.  The  new 
planning  tool  was  applied  to  study  a  hypothetical  phantom  case  and  a  prostate  case.  Compared  with  the 
conventional  inverse  planning  technique,  we  found  that,  for  the  same  target  dose  coverage,  the  critical 
structure  sparing  was  substantially  improved  for  all  testing  prostate  cases.  This  work  was  reported  as  an 
oral  presentation  in  2005  Annual  Meeting  of  American  Association  of  Physicists  in  Medicine  (AAPM)  at 
Seattle,  WA,  and  a  manuscript  reporting  the  details  of  the  work  is  in  progress. 

Evaluation  of  prostate  IMRT  dose  delivery  using  Kilovolt  (kV)  onboard  cone  beam  computed 
tomography  (CBCT):  kV  CBCT  based  on  flat-panel  technology  integrated  with  linear  accelerator  has 
recently  become  available  from  linac  vendors  for  therapy  guidance.  Currently,  the  system  is  primarily 
utilized  to  guide  the  patient  alignment.  As  an  advanced  tool  of  obtaining  a  patient’s  3D  representation, 
CBCT  also  affords  an  effective  means  for  us  to  examine  the  actual  dose  distribution  to  be  delivered  or 
already  delivered  to  the  patient  on  a  daily  basis.  We  have  recently  evaluated  the  accuracy  of  kV  CBCT- 
based  dose  calculation  and  addressed  some  logistic  issues  related  to  the  application  for  prostate  IMRT.  In 
reality,  image  quality  of  current  CBCT  is  not  as  good  as  conventional  diagnostic  CT  due  to  the  scatter  and 
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reconstruction  artifacts,  which  may  lead  to  significant  dosimetric  inaccuracy.  We  investigated  the 
feasibility  and  accuracy  of  CBCT-based  dose  calculation  and  to  proposed  a  deformable  electron  density 
mapping  (DEDM)  method  that  is  potentially  useful  to  facilitate  CBCT-based  dose  calculation.  In  the 
proposed  DEDM  technique,  the  CBCT  and  planning  CT  are  first  registered  by  using  a  deformable  image 
registration  model.  The  electron  density  distribution  is  then  mapped  from  the  planning  CT  to  the  CBCT. 
For  prostate  IMRT,  our  results  agree  with  the  planned  dose  distributions  to  within  3%.  The  use  of  DEDM 
reduces  the  dosimetric  inconsistency  between  the  CBCT-based  and  CT-based  dose  calculation  down  to  less 
than  1%  in  both  phantom  and  patient  studies.  Our  technique  provides  an  effective  way  to  ensure  that  the 
planned  IMRT  dose  distribution  can  be  realized  in  a  clinical  setting  and  should  have  significant  impact  on 
clinical  prostate  IMRT.  It  also  lay  the  foundation  for  future  adaptive  radiation  therapy.  This  work  has  been 
submitted  to  Medical  Physics  for  publication  (preprint  is  attached). 

Equivalent-uniform  dose  (EUD)-based  IMRT  Beam  Placement:  Beam  orientation  optimization  in 
IMRT  is  computationally  intensive  and  various  single  beam  ranking  techniques  have  been  proposed  to 
reduce  the  search  space.  Up  to  this  point,  none  of  the  existing  ranking  techniques  considers  the  clinically 
important  dose-volume  effects  of  the  involved  structures,  which  may  lead  to  clinically  irrelevant  angular 
ranking.  We  have  developed  a  clinically  sensible  angular  ranking  model  with  incorporation  of  dose- 
volume  effects  and  to  show  its  utility  for  IMRT  beam  placement.  The  general  consideration  in 
constructing  an  angular  ranking  function  is  that  a  beamlet/beam  is  more  preferable  if  it  can  deliver  a 
higher  dose  to  the  target  without  exceeding  the  tolerance  of  the  sensitive  structures  located  on  the  path  of 
the  beamlet/beam.  In  the  previously  proposed  dose-based  approach,  the  beamlets  are  treated 
independently  and,  to  compute  the  maximally  deliverable  dose  to  the  target  volume,  the  intensity  of  each 
beamlet  is  pushed  to  its  maximum  intensity  without  considering  the  values  of  other  beamlets.  When 
volumetric  structures  are  involved  (such  as  the  rectum  in  prostate  radiotherapy),  the  complication  arises 
from  the  fact  that  there  are  numerous  dose  distributions  corresponding  to  the  same  dose-volume  tolerance. 
In  this  situation,  the  beamlets  are  no  longer  independent  and  an  optimization  algorithm  is  required  to  find 
the  intensity  profile  that  delivers  the  maximum  target  dose  while  satisfying  the  volumetric  constraint(s).  In 
our  study,  the  behavior  of  a  volumetric  organ  was  modeled  by  using  the  EUD.  A  constrained  sequential 
quadratic  programming  algorithm  (CFSQP)  was  used  to  find  the  beam  profile  that  delivers  the  maximum 
dose  to  the  target  volume  without  violating  the  EUD  constraint(s).  It  is  shown  that  the  previously  reported 
dose-based  angular  ranking  represents  a  special  case  of  the  general  formalism  proposed  here.  We  also 
showed  that  the  proposed  technique  is  capable  of  producing  clinically  sensible  angular  ranking  for  a 
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variety  of  prostate  cases.  The  IMRT  plans  obtained  under  the  guidance  of  EUD-based  angular  ranking 
were  improved  in  comparison  with  that  obtained  using  the  conventional  uniformly  spaced  beams.  The 
EUD-based  function  is  a  general  approach  for  angular  ranking  and  allows  us  to  identify  the  potentially 
good  and  bad  angles  for  clinically  complicated  cases.  The  ranking  can  be  used  either  as  a  guidance  to 
facilitate  the  manual  beam  placement  or  as  prior  information  to  speed  up  the  computer  search  for  the 
optimal  beam  configuration.  Given  its  simplicity  and  robustness,  the  proposed  technique  should  have 
positive  clinical  impact  in  facilitating  the  IMRT  planning  process.  This  work  was  reported  as  an  oral 
presentation  in  2005  Annual  Meeting  of  American  Association  of  Physicists  in  Medicine  (AAPM)  at 
Seattle,  WA.  We  are  currently  refinining  the  technique  and  applying  it  to  generate  class-solution  for 
IMRT  prostate  treatment  planning. 

In  addition  to  the  above  three  projects,  we  are  extending  the  inverse  planning  infrastructure  for 
biologically  conformal  IMRT.  Currently,  an  inverse  planning  system  that  allows  us  to  utilize  the  spatial 
biology  distribution  does  not  exist.  A  theoretical  framework  to  quantitatively  incorporate  the  spatial 
biology  data  into  IMRT  inverse  planning  5  has  been  established.  Biologically  conformal  radiation  therapy 
(BCRT)  incorporating  patient  specific  biological  information  provides  an  outstanding  opportunity  for  us 
to  truly  individualize  radiation  treatment.  The  techniques  developed  in  this  proposal  will  find  natural 
application  in  the  next  generation  BCRT  treatment.  Finally,  the  Idea  Award  for  Prostate  Cancer  Research 
from  US  Amy  Medical  Research  and  Materiel  Command  also  provides  a  unique  educational  opportunity 
for  training  junior  researcher  through  the  participation  of  research  activities. 


III.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Established  a  theoretical  infrastructure  for  using  spatially  non-uniform  penalty  scheme  to  improve 
current  IMRT  inverse  planning. 

•  Developed  method  for  incorporating  existing  clinical  knowledge  into  inverse  planning  system. 

•  Proposed  an  electron  density  mapping  technique  based  on  deformable  image  registration  to  improve 
the  quality  of  cone-beam  CT  images. 

•  Established  a  robust  technique  for  using  onboard  cone-beam  CT  for  on-treatment  dose  validation  for 
prostate  IMRT. 

•  Improved  prostate  IMRT  beam  orientation  selection  procedure  by  integrating  organ  specific  EUD 
data. 
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•  Implemented  a  biological  model-based  BCRT  inverse  planning  algorithm  that  lays  the  technical 
foundation  for  next  generation  BCRT  treatment. 

IV.  REPORTABLE  OUTCOMES 

The  following  is  a  list  of  publications  resulted  from  the  grant  support  in  the  last  funding  period. 

Refereed  publications: 

1 .  Yang  Y,  Schreibmann  E.,  Li  T.,  King  C.,  and  Xing  L:  “Dosimetric  Evaluation  of  kV  Cone-Beam  CT- 
based  Dose  Calculation”.  Medical  Physics .  submitted,  2005. 

2.  Xing,  L,  Thorndyke  B,  Schreibmann  E,  Li  T,  Yang  Y,  Kim  G.,  Luxton  G,  Koong,  A,  Overview  of 
image  guided  radiation  therapy  (IGRT),  Medical  Dosimetry  (invited  review),  to  appear  in  March, 

2006. 

3.  Xia  P,  Yu  N,  Xing  L,  Syn  X,  Verhey,  L:  “Investigation  o  a  new  objective  function  for  inverse 
planning  optimization”.  Medical  Physics  32,  920-927,  2005. 

4.  Schreibmann  E  and  Xing  L:  “Image  registration  with  auto-mapped  control  volumes”.  Medical  Physics. 
Submitted. 

5.  Xing  L,  Levy,  D.  and  Yang  Y.,  Incorporating  clinical  outcome  data  into  inverse  treatment  planning, 
Medical  Physics,  manuscript  in  preparation. 

6.  Levy  D,  Paquin  D.,  Schreibmann  E.,  Xing  L,  Multiscale  registration  of  medical  imaging,  IEEE 
Transactions  on  Medical  Imaging,  to  be  submitted. 

Book  Chapter: 

1.  Xing  L.,  Yang  Y.,  Spielman  D.,  Molecular/Functional  Image-Guided  Radiation  Therapy,  T.  Bortfetld, 
R.  Schmidt-Ullrich,  W.  de  Neve  (editors),  Spinger-Verlag  Heidelberg,  Berlin. 

2.  Xing  L,  Wu  Q,  Yong  Y  and  Boyer  AL:  Physics  of  IMRT  and  Inverse  Treatment  Planning,  in  Intensity 
Modulated  Radiation  Therapy:  A  Clinical  Perspective,  Mundt  AF  and  Roeske  JC.  (editors),  page  20- 
51,  BC  Decker  Inc.  Publisher,  Hamilton,  Canada,  2005. 

3.  Wu  QW,  Xing  L,  Ezzel  G  and  Mohan  R,  Inverse  Treatment  Planning.”  In  Modern  Technology  of 
Radiation  Oncology.II.  Van  Dyk  J  (editor),  Medical  Physics  Publishing,  Madison,  WI,  2005. 

4.  Song  Y,  Boyer  A,  Ma  C,  Jiang  S.,  Xing  L,  Modulated  Electron  Beam  Therapy,  in  Intensity  Modulated 
Radiation  Therapy:  A  Clinical  Perspective,  Mundt  AF  and  Roeske  JC.  (editors),  page  327-336,  BC 
Decker  Inc.  Publisher,  Hamilton,  Canada,  2005. 

5.  Li  J  and  Xing  L,  Radiation  Dose  Planning,  Computer-Aided,  in  Encyclopedia  of  Medical  Devices  and 
Instrumentation,  John  G.  Webster  (editor),  John  Wiley  &  Sons,  in  press. 

Conference  abstract: 

1.  Xing,  L.  and  Spielman  D,  MRI/MRSI  and  Radiation  Therapy  Treatment  Planning,  2005  AAPM 
Annual  Meeting,  Seattle,  WA  (invited  talk). 

2.  Schreibmann  E.  and  Xing  L.,  EUD-based  beam  orientation  optimization,  oral  presentation  in  2005 
AAPM  Annual  Meeting,  Seattle,  WA. 

3.  Xing  L,  Levy,  D.  and  Yang  Y.,  Incorporating  clinical  outcome  data  into  inverse  treatment  planning, 
oral  presentation  in  2005  AAPM  Annual  Meeting,  Seattle,  WA. 

4.  Yang  Y,  Schreibmann  E.,  Li  T.,  King  C.,  and  Xing  L:  “Dosimetric  Evaluation  of  kV  Cone-Beam  CT- 
based  Dose  Calculation”,  oral  presentation  in  2005  AAPM  Annual  Meeting,  Seattle,  WA. 

5.  Yang  Y.  and  Xing  L,  Prescription  for  biologically  conformal  radiation  therapy,  poster  presentation  in 
2005  AAPM  Annual  Meeting,  Seattle,  WA. 
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6.  Schreibmann  E,  and  Xing  L:  “Image  registration  with  auto-mapped  control  volumes”.  Poster 
presentation  in  2005ASTRO  annual  meeting,  Denver. 

7.  Yang  Y  and  Xing  L:  “Optimization  of  radiation  dose -time-fractionation  scheme  with  consideration  of 
tumor  specific  biology”,  poster  presentation  in  2005ASTRO  annual  meeting,  Denver,  CO. 


V.  CONCLUSIONS 

An  infrastructure  has  been  established  to  execute  the  proposed  research.  Novel  IMRT  inverse 
planning  and  validation  techniques  are  being  developed  for  the  treatment  of  prostate  cancer.  A  few 
important  milestones  have  been  achieved  toward  the  general  goal  of  the  project.  These  include  (i) 
established  a  theoretical  infrastructure  of  spatially  non-uniform  penalty  scheme  for  inverse  planning;  (ii) 
developed  method  for  incorporating  existing  clinical  knowledge  into  inverse  planning;  (iii)  proposed  an 
electron  density  mapping  technique  to  improve  the  quality  of  cone-beam  CT  (CBCT)  images;  (iv) 
established  a  robust  technique  for  using  onboard  CBCT  for  on-treatment  IMRT  dose  validation;  (v) 
improved  prostate  IMRT  beam  orientation  selection  by  integrating  organ  specific  EUD.  We  have  also 
implemented  a  biological  model-based  BCRT  inverse  planning  algorithm  that  lays  the  technical  foundation 
for  next  generation  BCRT  treatment.  Integration  and  further  refinement  of  the  above  tools  are  underway.  It 
is  expected  these  tools  will  substantially  improve  the  current  prostate  IMRT  treatment. 
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Abstract 

Kilovolt  (kV)  CBCT  based  on  flat-panel  technology  integrated  with  linear  accelerator  has 
recently  become  available  from  linac  vendors  for  therapy  guidance.  Currently,  the  system 
is  primarily  utilized  to  guide  the  patient  alignment.  As  an  advanced  tool  of  obtaining  a 
patient’s  3D  representation,  CBCT  also  affords  an  effective  means  for  us  to  examine  the 
actual  dose  distribution  to  be  delivered  or  already  delivered  to  the  patient  on  a  daily  basis. 
Before  this  can  be  implemented  clinically,  the  accuracy  of  kV  CBCT-based  dose 
calculation  must  be  evaluated  and  some  logistic  issues  related  to  the  application  need  to 
be  addressed.  Indeed,  image  quality  of  current  CBCT  is  not  as  good  as  conventional 
diagnostic  due  to  the  scatter  and  organ  motion  artifacts,  which  lead  to  significant 
dosimetric  inaccuracy.  This  work  is  aimed  to  investigate  the  feasibility  and  accuracy  of 
CBCT-based  dose  calculation  and  to  propose  a  deformable  electron  density  mapping 
(DEDM)  method  that  is  potentially  useful  to  facilitate  CBCT-based  dose  calculation.  In 
the  proposed  DEDM  technique,  the  CBCT  and  planning  CT  are  first  registered  by  using  a 
deformable  image  registration  model.  The  electron  density  distribution  is  then  mapped 
from  the  planning  CT  to  the  CBCT.  The  CBCT  with  .the  mapped  electron  density 
information  serves  as  the  backbone  for  more  accurate  CBCT-based  dose  calculation.  For 
disease  sites  where  intra-fractional  organ  motion  is  not  an  issue,  this  study  indicates  that 
CBCT  can  be  employed  directly  for  dose  calculation  without  relying  on  the  DEDM  and 
the  results  agree  with  the  planned  dose  distributions  to  within  3%.  On  the  other  hand,  in 
the  presence  of  motion  artifacts,  the  Hounsfield  number  distribution  can  altered 
substantially  and,  as  a  result,  it  is  found  that  the  dosimetric  inaccuracy  can  be  more  than 
5%  when  the  bare  CBCT  is  used.  The  use  of  DEDM  reduces  the  dosimetric  inconsistency 
between  the  CBCT-based  and  CT-based  dose  calculation  down  to  less  than  1%  in  both 
phantom  and  patient  studies.  While  the  true  solution  to  the  hurdle  lies  in  the  effective 
removal  of  motion  artifacts  in  CBCT,  the  DEDM  approach  seems  to  afford  a  useful 
interim  technique  for  improved  CBCT-based  dose  calculation. 

Key  word:  CBCT,  IGRT,  Dose  verification.  Deformable  registration 
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I.  Introduction 

Modern  radiation  therapy  techniques,  such  as  3D  conformal  radiotherapy  (3DCRT)  and 
intensity-modulated  radiation  therapy  (IMRT),  provide  unprecedented  means  for 
producing  exquisitely  shaped  radiation  doses  that  closely  conform  to  the  tumor 
dimensions  while  sparing  sensitive  structures.  As  a  result  of  greatly  enhanced  dose 
conformality,  more  accurate  beam  targeting  becomes  an  urgent  issue  in  radiation  therapy. 
In  current  practice,  large  uncertainties  exist  in  tumor  target  localization  due  to  intra-  and 
inter-organ  motions  during  the  course  of  radiation  treatment.  As  thus,  large  safety 
margins  around  the  tumor  targets  and  sensitive  structures  are  introduced  to  cope  with  the 
otherwise  insoluble  patient  localization  problem.  The  use  of  non-optimal  margins 
compromises  the  patient  care  and  adversely  affects  the  treatment  outcome  ^  "6.  The  need 
to  improve  targeting  in  high  precision  radiation  therapy  has  recently  spurred  a  flood  of 
research  activities  in  image-guided  radiation  therapy  (IGRT)7_1 1. 

CBCT  based  upon  flat-panel  technology  integrated  with  a  medical  linear 
accelerator  has  recently  become  available  from  Linac  vendors  for  therapy  guidance.  The 
volumetric  images  may  be  used  to  verify  and  correct  the  planning  patient  setup  in  the 
linac  coordinates  by  comparing  with  the  patient  position  defined  in  treatment  plan.  Both 

kV  and  MV  beams  *  2,  13  have  been  utilized  for  this  application.  The  former  typically 
consists  of  a  kV-source  and  flat-panel  combination  mounted  on  the  drum  of  a  medical 

accelerator  with  the  kV  imaging  axis  orthogonal  to  that  of  MV  therapy  beam.  In 
addition  to  guide  the  patient  setup  process,  CBCT  data  acquired  prior  to  the  treatment 
can,  in  principle,  be  used  to  recalculate  or  verify  the  treatment  plan  for  the  patient 
anatomy  of  the  treatment  day.  The  recalculation  starts  with  the  intended  fluence  maps 
from  the  patient’s  treatment  plan,  whereas  the  verification  is  done  by  using  the  fluence 
maps  measured  at  the  exiting  sides  of  the  incident  beams.  If  CBCT-based  dose 
calculation  is  accurate  enough  (say,  with  an  accuracy  within  3%),  this  will  provide  a 
valuable  option  for  us  to  predict/assess  the  patient  dose  on  a  daily  basis.  In  reality, 
because  of  the  presence  of  organ  movement/deformation,  it  is  conceivable  that  the  dose 
distributions  delivered  to  the  patient  are  usually  different  from  fraction  to  fraction.  It  is 
paramount  to  be  able  to  monitor  the  actual  patient  dose  for  each  fraction  as  well  as  the 
accumulated  doses  to  the  target  and  sensitive  structures  while  the  fractional  treatment 
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proceeding.  This  will  not  only  give  physician  more  confidence  about  the  treatment  but 
may,  in  future,  afford  us  an  effective  means  to  adaptively  modify  the  patient’s  treatment 
plan  during  the  course  of  a  radiation  therapy  based  on  the  dose  that  has  already  been 
delivered. 

The  accuracy  of  MV  fan-beam  and  cone-beam  CT  has  recently  been  assessed  by 
Langen  et  a\^  and  Poulliot  et  al  The  potential  of  its  counterpart,  the  kV  CBCT,  for 
dosimetric  calculation  has,  on  the  other  hand,  not  been  examined  systematically.  The 
purpose  of  this  work  is  two-fold:  to  evaluate  the  dosimetric  accuracy  of  kV  CBCT-based 
dose  calculation  and  to  explore  a  deformable  electron  density  mapping  strategy  for 
improving  the  performance  of  the  calculation.  For  disease  sites  where  intra-fractional 
organ  motion  is  not  an  issue,  the  general  reference  drawn  from  this  study  is  that  it  is 
adequate  to  directly  use  current  CBCT  for  dose  calculation  because  the  dosimetric 
inaccuracy  is  generally  less  than  3%  as  compared  with  that  based  upon  the  conventional 
CT.  On  the  other  hand,  when  dealing  with  regions  in  the  thorax  or  upper  abdomen,  the 
motion  artifacts  in  CBCT  are  found  quite  severe  and  significant  dosimetric  errors  (~5%) 
are  observed.  In  this  situation,  the  deformable  electron  density  mapping  method  seems  to 
be  valuable  for  more  accurate  CBCT-based  dose  recalculation. 


II.  Method  and  Materials 

A.  Data  acquisition 

The  on-board  imager  (OBI)  integrated  in  a  Trilogy™  medical  linear  accelerator  (Varian 
Medical  Systems,  Palo  Alto,  CA)  is  used  in  this  work  to  acquire  CBCT  images.  The  kV 
OBI  system  is  capable  of  obtaining  low-dose,  high-resolution  radiography,  fluoroscopy 
and  CBCT.  The  system  is  mounted  on  the  treatment  machine  via  robotically  controlled 
arms,  which  operate  along  three  axes  of  motion.  It  can  be  automatically  positioned  using 
the  robotic  technology  and  control  software,  making  OBI  clinically  very  practical.  A  150 
kV  X-ray  tube  with  maximum  32  ms  pulse  length  for  continuous  irradiation  and 
maximum  320  ms  pulse  length  for  single  pulse  is  designed  for  generating  high-resolution 
images  from  a  moving  gantry.  The  spot  of  the  tube  is  located  at  90°  to  the  MV  source  and 
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100cm  from  the  radiation  axis  of  the  accelerator.  A  39.7cm  X  29.8  cm  amorphous  silicon 
flat-panel  X-ray  image  detector  (Varian  PortalVision™  aSlOOO)  mounted  opposite  the 
kV  tube  is  used  to  acquire  digital  images  with  a  pixel  matrix  of  2048  X  1536.  Using  the 
OBI  system,  the  CBCT  data  can  be  acquired  in  two  modes:  full  fan  mode  and  half  fan 
mode.  In  the  full  fan  mode,  675  projections  are  taken  during  the  whole  364°  gantry 
rotation  with  a  field  of  view  (FOV)  about  26.6  cm  in  diameter  and  17cm  in  length.  The 
data  acquisition  time  is  about  45  second  and  the  reconstruction  time  for  340  slices  of 
512X512  CBCT  images  with  a  voxel  size  of  0.5mm3  is  also  about  1  minute  on  a  PC.  The 
half  fan  mode  is  designed  to  obtain  larger  FOV.  In  this  mode,  about  965  projections  are 
taken  during  the  364°  gantry  rotation  and  a  FOV  of  45  cm  diameter  X  15cm  can  be 
achieved.  The  data  acquisition  and  reconstruction  time  for  512X512  CBCT  images  with  a 
voxel  size  of  0.95mm’  using  this  mode  is  about  double  compared  with  the  full  fan  mode. 

B.  Calibration  of  relative  electron  density 

To  use  CT  or  CBCT  for  radiation  dose  calculation,  it  is  required  to  relate  the  Flounsfield 
Unit  (HU)  of  the  scanner  with  the  actual  electron  density.  A  CT-phantom,  Catphan-600 
module  CTP404  (Phantom  Laboratory,  NY),  was  used  for  the  calibration  of  planning  CT 
(GE  Discovery-ST  PET/CT  scanner,  Milwaukee,  WI)  and  CBCT.  The  CTP404  has  a 
diameter  of  150  mm  and  contains  17  different  sizes  of  inserts  with  seven  different  tissue 
substitute  materials,  air,  PMP,  LDPE,  Polystryrene,  Acrylic,  Delrin  and  Teflon, 
respectively.  Their  relative  electron  densities  ranged  from  0  to  1 .866.  A  cross  section  of 
the  phantom  is  shown  in  figure  2.  The  calibration  of  a  CT  scanner  involves  acquiring  CT 
images,  obtaining  average  HUs  for  each  inserting  materials,  and  plot  the  HU  as  a  function 
of  the  relative  electron  density.  For  CBCT,  it  is  necessary  to  calibrate  separately  for  full 
and  half  fan  modes  because  the  beam  geometry  and  characteristics  of  the  two  types  of 
scanning  modes  are  different.  More  details  on  the  calibration  procedure  can  be  found  in 
the  manual  of  Trilogy. 

In  order  to  test  the  stability  of  the  CBCT  calibration  curve  with  time,  the  phantom 
was  repeatedly  scanned  every  week  for  two  months.  The  obtained  HU  vs  relative  electron 
density  curves  were  compared  to  assess  the  HU  fluctuations  with  time. 
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C.  Phantom  study 

CT  and  CBCT  images  of  the  Catphan-600  phantom  were  acquired  using  the  procedure 
outlined  in  Sec.  II.A.  The  phantom  was  placed  on  a  platform  that  can  be  set  to  one¬ 
dimensional  cyclic  motion  with  a  number  of  speeds.  The  axis  of  the  cylindrically  shaped 
phantom,  along  which  the  phantom  moves  cyclically,  was  angled  from  the  central  axis  of 
the  couch  plane  by  about  15°  in  order  for  the  on-board  imager  to  “see”  or  “realize”  the 
motion  when  the  cylindrical  phantom  moves.  The  movement  of  the  phantom  produces 
motion  artifacts  in  the  images  and  allows  us  to  evaluate  the  performance  of  CBCT-based 
dose  calculation  in  the  presence  of  organ  motion.  The  full  fan  mode  was  used  to  scan  the 
phantom.  CT  and  CBCT  images  of  the  phantom  were  acquired  with  and  without  motion. 
In  the  former  case,  the  peak-to-peak  amplitude  of  the  motion  was  0.5  cm  and  the  period 
was  4s,  which  mimics  the  situation  of  a  patient’s  breathing  motion. 

To  quantify  the  difference  in  the  image  quality  of  the  CT  and  CBCT  images,  we 
first  analyzed  the  HU  distribution  for  the  four  sets  of  CT  images,  corresponding  to 
planning  CT  and  CBCT  with  and  without  motion.  The  influence  of  phantom  motion  on 
the  HU  distribution  was  investigated.  The  CT  and  CBCT  images  were  imported  to  a 
Varian  Eclipse  treatment  planning  system  for  dosimetric  comparison  study.  For  planning 
and  evaluation  purpose,  a  hypothetical  spherical  target  with  a  diameter  of  5cm  was 
created  at  the  center  of  the  phantom  and  a  single  5  X  5cm2  6MV  photon  beam  was  used 
to  irradiate  the  target.  A  simple  beam  configuration  was  used  here  because,  in  this  way, 
the  results  are  more  intuitively  interpretable.  Four  plans,  corresponding  to  the  four 
different  sets  of  CT  images,  were  generated  using  the  same  target  and  beam 
configuration.  The  pencil  beam  convolution  dose  calculation  algorithm  implemented  in 
Varian  Eclipse  treatment  planning  system  was  adopted  for  dose  calculation.  The  resultant 
isodose  curves,  dose  profiles  and  DVHs  were  compared. 

D.  Patient  study 

A  prostate  cancer  patient  and  a  lung  cancer  patient  were  selected  for  the  evaluation  study 
of  CBCT-based  dose  calculation  and  to  demonstrate  the  feasibility  of  the  proposed 
deformable  electron  density  mapping  (DEDM)  technique  (see  next  sub-section)  for 
improved  dose  calculation  accuracy.  For  the  prostate  case,  the  targets  included  the  PTV, 
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consisting  of  the  prostate  gland  with  a  margin  of  6mm  and  the  seminal  vesicles.  The 
critical  structures  were  rectum,  bladder  and  femoral  heads.  An  IMRT  plan  using  five 
15MV  photon  beams  with  gantry  angles  of  35°,  110°,  180°,  250°,  and  325°  (in  1EC 
convention)  was  adopted  for  the  prostate  case.  The  plan  was  normalized  to  deliver  a 
prescription  dose  of  78Gy  to  99%  the  prostate  PTV  and  no  less  than  50Gy  to  the  98%  of 
seminal  vesicles  in  39  fractions.  For  the  lung  cancer  case,  the  PTV  consisted  of  the  CTV 
with  a  margin  of  10mm.  The  critical  structures  involved  were  the  right  lung  and  the 
spinal  cord.  A  conventional  3D  conformal  plan  with  three  15MV  photon  beams  (45°, 
180°  and  288°  in  IEC  convention)  was  generated  for  this  case.  In  this  plan,  the  field  shape 
of  each  beam  was  determined  by  conforming  the  PTV  projection  in  the  corresponding 
beam  direction.  The  plan  was  normalized  to  deliver  a  prescription  dose  of  70Gy  to  100% 
of  the  target  volume  in  35  fractions.  The  CBCT  images  of  the  patients  were  obtained 
using  the  half  fan  mode.  The  CT  and  CBCT  images  were  registered  using  a  rigid  image 
registration  package  provided  in  the  Eclipse  treatment  planning  system.  For  each  case, 
the  IMRT  or  3D  CRT  planning  parameters,  which  include  beam  configuration,  MU 
settings,  and  MLC  files,  used  for  treatment  were  employed  to  recalculate  the  dose  based 
on  the  CBCT  data.  The  CT  and  CBCT-based  treatment  plans  were  then  compared. 

E.  Deformable  electron  density  mapping 

The  dosimetric  inaccuracy  of  CBCT-based  dose  calculation  in  thorax  and  abdomen  arises 
from  the  inability  of  the  CBCT  technique  to  provide  accurate  FIU  or  relative  electron 
density  distribution,  primarily  due  to  organ  motion  induced  artifacts.  The  genuine 
solution  to  the  problem  lies  in  the  improvement  of  the  CBCT  acquisition  technology  so 
that  artifact-free  images  can  be  acquired.  While  this  endeavor  is  still  on-going,  here  we 
propose  an  interim  solution  for  dealing  with  the  problem.  Under  the  assumptions  that  the 
HU  or  relative  electron  density  distribution  is  known  from  planning  CT  and  an  acceptable 
geometric  registration  between  CT  and  CBCT  is  achievable  by  a  deformable  registration, 
we  propose  to  map  the  electron  density  in  the  planning  CT  onto  the  daily  setup  CBCT 
and  then  carry  out  the  dose  calculation.  The  CBCT  with  mapped  electron  density, 
referred  to  as  modified  CBCT,  possesses  the  anatomical  information  of  CBCT  and  yet  the 
electron  density  information  of  the  planning  CT.  Dose  calculation  based  on  the  modified 
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CBCT  allows  us  to  compute  more  accurately  the  delivered  dose  with  the  patient  in  his/her 
setup  position.  The  mapping  process  is  described  as  follows. 

A  free  form  spline  (BSpline)  deformable  model  15,  16  was  employed  to  register 
the  planning  CT  and  CBCT  and  map  the  deformed  electron  density  from  planning  CT  to 
CBCT.  The  method  was  used  for  several  IGRT  related  projects  in  our  group  and  others 
and  its  simplicity  and  accuracy  have  been  demonstrated  1^>  18.  Briefly,  in  the  BSpline 
model,  a  lattice  of  user-defined  nodes  is  overlaid  on  the  image.  Each  node  contains  a 
deformation  vector,  whose  components  are  to  be  determined  by  optimizing  a  metric 
function  that  characterizes  the  goodness  of  the  registration.  The  variables  of  the  metric 
function  consisted  of  the  coordinates  of  the  BSpline  nodes.  In  this  work,  a  voxel-based 
normal  cross  correlation  (NCC)  metric  was  used.  A  suitable  set  of  node  deformations  was 
determined  using  the  gradient-based  algorithm  L-BFGS  16,  19;  which  is  known  for  its 
superior  performance  in  large-scale  optimization  problems.  The  optimizer  iteratively 
varies  the  nodal  displacements  to  optimize  the  metric.  The  deformation  at  any  point  of 
the  image  is  calculated  by  spline  interpolation  of  closest  nodes  values.  Unlike  other  spline 
models,  the  BSplines  are  locally  controlled,  such  that  the  displacement  of  an  interpolation 
point  is  influenced  only  by  the  closest  grid  points  and  changing  a  lattice  node  only  affects 
the  transformation  regionally,  making  it  efficient  in  describing  local  deformations.  After 
the  deformable  registration,  the  HU  in  each  voxel  in  planning  CT  was  mapped  to  the 
corresponding  point  in  the  reference  CBCT  to  produce  the  modified  CBCT  images. 

The  feasibility  of  DEDM  technique  was  evaluated  by  using  the  two  patients 
mentioned  above.  For  this  purpose,  the  CT  and  CBCT  images  were  registered  using  the 
BSpline  model.  The  targets  and  sensitive  structures  contoured  on  the  planning  CT  were 
copied  to  the  CBCT  using  the  deforfnable  model.  For  each  patient,  the  treatment  plan 
parameters  were  employed  to  recompute  the  dose  distribution  based  on  the  patient’s 
modified  CBCT.  The  resultant  isodose  curves  and  DVHs  were  evaluated  and  the  level  of 
improvement  in  dosimetry  due  to  the  use  of  DEDM  was  assessed. 


III.  Results 
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A.  Calibration  of  CT  and  CBCT 

The  relation  between  kV  HU  distribution  of  CBCT  and  relative  electron  dosimetry  was 
established  by  using  a  Catphan-600  CT  phantom  following  the  procedure  described  in 
Sec.  II.B.  The  calibration  curves  for  planning  CT,  half  fan  and  full  fan  CBCT  modes  are 
shown  in  figure  la.  Figure  lb  compares  the  calibration  curves  obtained  with  an  interval 
of  1  week  during  continuous  two  months  for  full  fan  CBCT.  No  significant  fluctuations 
were  found  in  the  calibration  data,  which  is  similar  to  what  have  been  observed  for 

MV  14  stability  of  the  kV  CBCT  and  electron  density  calibration  is  good  indicator  of 
the  HU  number  integrity  and  the  overall  performance  of  the  CBCT  system. 

B.  Phantom  study 

Figures  2a  to  2d  show  the  same  transverse  slices  of  the  CT  and  CBCT  images  of  the 
Catphan-600  phantom  with  and  without  motion.  The  first  two  panels  are  the  CT  and 
CBCT  images  of  the  phantom  in  the  absence  of  motion,  and  the  second  two  show  the 
same  with  the  phantom  motion  “switched  on”.  It  is  seen  that  the  quality  of  CBCT  images 
is  worse  than  that  of  the  conventional  planning  CT,  especially  in  the  presence  of  motion. 
The  HU  profiles  of  the  four  images  along  the  two  orthogonal  lines  (lines  A-A  and  B-B  as 
marked  in  figure  2)  are  plotted  in  figure  3.  It  is  found  that  the  HU  profiles  of  the  planning 
CT  and  CBCT  normally  agree  to  within  10%  in  the  static  situation.  On  the  other  hand, 
when  the  motion  is  “switched  on”,  CBCT  shows  a  much  greater  level  of  artifacts  (figure 
2d)  and  the  HU  difference  between  the  conventional  CT  and  CBCT  is  aggravated,  with 
the  maximum  difference  reaching  several  hundred  HUs. 

Figures  4,  5  and  6  present  the  dosimetric  results  calculated  using  a  single  6  MV  5 
X  5cm2  photon  beam.  Figures  4a  to  4d  depict  the  dose  distributions  in  a  transverse  slice 
calculated  based  on  the  four  sets  of  images  given  in  figure  2.  Figures  5a  and  5b  compare 
the  dose  profiles  along  the  two  orthogonal  lines  (lines  A-A  and  B-B  in  figure  2),  and 
figure  6  compares  the  DVHs  of  the  target  for  the  four  different  situations.  From  these 
results  we  find  that  the  dose  calculated  using  planning  CT  agrees  with  that  of  CBCT- 
based  calculation  to  within  1 .0%,  indicating  that  it  is  acceptable  to  use  kV  CBCT  for  dose 
calculation  if  no  organ  motion  presents.  However,  when  phantom  motion  is  involved,  the 
motion  induced  artifacts  significantly  influence  the  HU  distribution  and  thus  the  accuracy 
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of  CBCT-based  dose  calculation.  For  this  simple  phantom  case,  we  find  that  the 
discrepancy  between  the  planning  CT-  and  CBCT-based  calculations  is  about  3%,  which 
is  clinically  significant.  The  motion  artifacts  existing  in  current  CBCT  limit  the  direct  use 
of  CBCT  for  dose  calculation  when  intra-fractional  organ  motion  is  not  negligible. 

C.  Patient  study 

Figures  7a  to  7c  show  the  same  transverse  slices  of  the  planning  CT,  CBCT,  and 
checkerboard  image  resulting  from  the  deformable  registration  of  the  two  sets  of  images 
for  the  prostate  case.  The  modified  CBCT  obtained  by  mapping  the  HUs  from  the 
planning  CT  to  CBCT  is  shown  in  figure  7d.  Our  previous  studies  have  indicated  that  a 
registration  accuracy  better  than  2mm  is  achievable  by  using  the  BSpline  deformable 

registration  technique^  described  in  the  last  section.  As  can  be  seen  from  the 
checkerboard  overlay,  a  good  registration  between  CT  and  CBCT  is  obtained.  Figure  8 
shows  the  isodose  distributions  for  the  three  types  of  calculations  based  on  planning  CT, 
CBCT,  and  modified  CBCT.  A  comparison  of  DVHs  of  PTV,  prostate,  seminal  vesicles, 
bladder  and  rectum  is  presented  in  figure  9.  The  results  demonstrate  that,  while  there  is 
significant  dosimetric  discrepancy  between  the  results  obtained  based  on  the  planning  CT 
and  modified  CBCT,  the  results  obtained  using  the  CBCT  and  modified  CBCT  is  similar, 
except  for  the  seminal  vesicles,  in  which  the  DVH  difference  is  somewhat  large  due  to  its 
relatively  small  volume  and  more  deformation.  Similar  phenomenon  was  also  observed 
by  the  MD  Anderson  group  using  their  daily  CT  on-rail^,  20  [n  general,  the  difference 
between  the  planned  dose  distribution  and  that  computed  based  on  CBCT  arises  from  two 
factors:  (i)  patient  positioning  error  and  organ  deformation/displacement;  and  (ii)  relative 
electron  density  difference  between  the  two  sets  of  CT  images.  The  small  discrepancy 
between  the  doses  computed  using  CBCT  and  modified  CBCT  suggests  that,  in  the 
prostate  case,  the  second  factor  is  small  and  it  is  acceptable  to  directly  use  CBCT  for  dose 
calculation.  The  dosimetry  is  predominantly  determined  by  the  accuracy  of  patient  setup 
and  the  level  of  interfractional  deformation  /displacement  of  the  prostate,  rectum  and 
bladder. 

Figure  10a  to  lOd  show  the  same  transverse  slice  from  the  planning  CT,  CBCT, 
checkerboard  overlay  of  the  planning  CT  and  the  modified  CBCT  for  the  lung  cancer 
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case.  Figures  1 1  and  12  compare  the  isodose  distributions  and  the  DVHs  of  the  target  and 
sensitive  structures  calculated  based  on  the  three  sets  of  images  for  the  three-beam  3D 
conformal  treatment  plan.  In  this  case,  the  dosimetric  discrepancy  between  the  CBCT- 
and  modified  CBCT-based  calculations  is  much  larger  than  that  in  the  prostate  case, 
especially  for  the  right  lung  and  PTV.  The  maximum  dose  difference  is  about  5%. 
However,  the  discrepancy  between  the  results  obtained  using  planning  CT  and  modified 
CBCT  becomes  much  less.  This  study  exemplifies  that,  unlike  the  prostate  case,  the 
dosimetric  inaccuracy  arising  from  the  inferior  image  quality  of  CBCT  in  the  dose 
verification  calculation.  The  motion  artifacts  not  only  make  it  difficult  to  see  the  extent  of 
the  tumor,  but  also  limit  the  direct  use  of  CBCT  for  dose  calculation. 


IV.  Discussion 

The  feasibility  and  accuracy  of  using  kV  CBCT  to  calculate  dose  have  been  investigated 
with  a  phantom  and  two  clinical  cases.  In  the  absence  of  motion  artifacts,  it  seems  to  be 
acceptable  to  directly  use  CBCT  for  dose  verification  calculation.  Otherwise,  extra 
caution  is  required  to  avoid  significant  dosimetric  inaccuracy.  To  cope  with  the  problem 
of  deteriorated  imaging  quality  of  CBCT  in  the  presence  of  organ  motion,  a  DEDM 
method  has  been  proposed  to  map  the  electron  density  information  from  the  patient’s 
planning  CT  to  the  setup  CBCT  with  a  deformable  image  registration.  In  IGRT,  since  the 
registration  has  to  be  done  for  the  purpose  of  patient  setup,  the  computational  overhead  of 
introducing  DEDM  is  minimal.  Before  an  effective  CBCT  motion  artifacts  removal 
technique  is  in  place,  DEDM  provides  a  useful  interim  solution  to  the  problem. 

Dose  distributions  computed  based  on  CBCT  and  modified  CBCT  represent  the 
dose  to  be  delivered  to  the  patient  because  the  CBCT  was  acquired  prior  to  the  patients’ 
treatments  after  the  patients  were  repositioned/shifted  using  the  patient  setup  procedure  in 
current  practice.  In  the  prostate  IMRT  plan,  the  inherent  dosimetric  error  resulted  from 
the  use  of  CBCT  images  is  small.  However,  the  dosimetric  error  caused  by  the  inter- 
fractional  organ  motion/deformation  is  not  insignificant,  as  revealed  by  using  dose  re¬ 
calculation  results  given  in  the  last  section.  A  few  groups  are  working  on  deformable 
model  based  segmentation  and  patient  setup  procedures^,  9,  21,  22  when  deformable 
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registration  is  used,  there  are  a  few  options  to  achieve  the  registration  depending  on 
whether  the  primary  aim  is  to  match  soft-tissue,  or  to  align  3D  bony  structures.  The 
multiple  choices  resulting  from  the  fact  that  the  dimensionality  of  the  patient  data  is  much 
greater  than  that  in  the  patient  setup  procedure  and  suggest  that  deformable  registration  is 
not  the  ultimate  solution  to  volumetric  image-guided  radiation  therapy.  The  technique 
improves  the  current  body-structure-based  patient  alignment  method  since  it  partially 
takes  into  account  organ  deformation  by  achieving  the  closest  overlay  match  possible 
between  the  planning  and  CBCT  data  sets  according  to  our  clinical  objective,  and  provide 
an  improved  potion  positioning  technique.  We  should  emphasize  that,  even  when  3D 
volumetric  based  deformable  registration  is  available  in  the  future,  the  problem  of  patient 
positioning  will  not  disappear  as  relative  organ  deformations/displacements  may  well 
persist.  A  possible  solution  to  accommodate  various  factors  mentioned  above  is  for  us  to 
re-optimize  or  tweak  the  IMRT  plan  based  on  the  patient’s  setup  CBCT.  Indeed,  in  order 
to  fully  utilize  the  CBCT  volumetric  data,  a  new  paradigm  with  seamlessly  integrated 
simulation,  planning,  verification,  and  delivery  procedure  is  urgently  needed.  Until  this  is 
realized  clinically,  the  volumetric  imaging  is  nothing  but  an  expensive  extension  of  the 
existing  planar  verification  approach. 

For  the  lung  3D  conformal  plan,  we  found  that  the  dosimetric  error  for  the  target 
from  the  organ  interfractional  motion/deformation  and  intrafractional  motion  is  small. 
The  main  reasons  are  that  the  organ  interfractional  motion/deformation  is  not  as 
significant  as  that  in  the  prostate  case  and  enough  margins  are  given  in  the  3D  conformal 
plan.  Furthermore,  the  dose  distribution  in  the  target  is  more  uniform  in  this  case. 
However,  the  dosimetric  error  from  the  lower  quality  of  the  CBCT  images  becomes 
larger  since  the  influence  of  the  intrafractional  motion  on  the  CBCT  image  quality.  In  this 
situation,  it  is  suggested  that,  in  order  to  accurately  evaluate  the  delivered  dose 
distributions,  the  deformed  CT  images  generated  using  the  proposed  deformable  electron 
density  mapping  technique  are  required  to  replace  the  original  CBCT  images  for  the  task 
of  dose  verification  calculation. 


V.  Conclusion 
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On-board  CBCT  provides  useful  volumetric  anatomy  information  for  patient  positioning 
verification.  When  used  for  dose  verification  calculation,  it  is  required  to  have  a  reliable 
HU  to  electron  density  curve.  Our  phantom  and  patient  study  have  indicated  that,  in  the 
absence  of  motion  artifacts,  the  dosimetric  inaccuracy  is  generally  less  than  1%  as 
compared  with  the  calculation  of  the  treatment  planning  system.  The  dosimetric  errors 
are  much  more  pronounced  when  intra-fractional  organ  motion  is  present.  In  this 
situation,  a  direct  use  of  CBCT  for  dose  calculation  is  not  recommended.  The  use  of  a 
deformable  registration  permits  us  to  incorporate  the  electron  density  distribution  from 
the  planning  CT  and  to  calculate  the  dose  more  accurately.  The  proposed  DEDM 
approach  affords  a  practical  solution  to  estimate  the  dose  to  be  delivered  or  already 
delivered  to  the  patient  based  on  the  setup  CBCT. 
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Figure  Captions 


Fig.  1.  (a)The  calibration  curves  (Hounsfield  number  vs  relative  electron  density)  for 
planning  CT,  half  fan  and  full  fan  mode  CBCT;  (b)  the  variation  of  calibration  curves 
with  time  for  the  full  fan  CBCT. 

Fig.  2.  The  CT  and  CBCT  images  with  and  without  motion  for  the  Catphan-600:  (a) 
planning  CT  in  the  absence  of  phantom  motion;  (b)  CBCT  in  the  absence  of  phantom 
motion;  (c)  planning  CT  with  moving  phantom;  and  (d)  CBCT  with  moving  phantom. 

Fig.  3.  HU  profiles  of  planning  CT  and  CBCT  images  (see  figure  2)  along  the  A-A  line 
(panel  a)  and  B-B  line  (panel  b). 

Fig.  4.  Dose  distributions  in  a  transverse  slice  calculated  based  on  the  four  sets  of  CT  data 
shown  in  figure  2:  (a)  planning  CT;  (b)  CBCT;  (c)  planning  CT  with  a  motion;  and  (d) 
CBCT  with  a  motion.  In  all  four  situations,  a  5  X  5cm2  single  field  plan  was  used  to 
irradiate  a  spherical  hypothetical  target  with  a  diameter  of  5cm  located  at  the  phantom 
center. 

Fig.  5.  Comparison  of  the  dose  profiles  along  the  two  orthogonal  lines  shown  in  figure  2 
for  the  Catphan-600  phantom:  (a)  profile  along  the  A-A  line;  (a)  profile  along  the  B-B 
line. 

Fig.  6.  Comparison  of  the  target  DVHs  calculated  based  on  the  four  sets  of  CT  data 
shown  in  figure  2  for  the  phantom  case. 

Fig.  7.  CT,  CBCT  and  modified  CBCT  images  for  the  prostate  case:  (a)  planning  CT;  (b) 
daily  CBCT;  (c)  checkerboard  overlay  of  CT  and  CBCT  after  the  deformation 
registration;  and  (d)  modified  CT. 
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Fig.  8.  Dose  distributions  in  a  transverse  slice  calculated  based  on  the:  (a)  planning  CT; 
(b)  CBCT;  and  (c)  modified  CB  CT  for  the  prostate  case. 

Fig.  9.  DVHs  of  the  prostate,  PTV,  rectum  and  bladder  obtained  based  on  the  planning 
CT,  CBCT  and  modified  CBCT  images  for  the  prostate  case. 

Fig.  10.  CT,  CBCT  and  modified  CT  images  for  the  lung  case:  (a)  planning  CT;  (b) 
CBCT;  (c)  checkerboard  image  after  the  deformation  registration;  and  (d)  modified 
CBCT. 

Fig.  1 1.  Dose  distribution  in  a  transverse  slice  calculated  based  on:  (a)  planning  CT;  (b) 
CBCT;  and  (c)  modified  CBCT  for  the  lung  case. 

Fig.  12.  Comparison  of  the  DVFls  of  GTV,  PTV,  right  lung  and  spinal  cord  obtained 
based  on  planning  CT,  CBCT  and  the  modified  CBCT  images. 
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Abstract 

A  multiscale  image  registration  technique  is  presented  for  the  registration  of  medical  images  that 
contain  significant  levels  of  noise.  Registration  is  achieved  by  obtaining  a  hierarchical  multiscale 
decomposition  of  the  noisy  images  and  registering  the  resulting  components.  This  approach  enables 
successful  registration  of  images  that  contain  noise  levels  well  beyond  the  amount  at  which  ordi¬ 
nary  registration  fails.  Experiments  are  presented  that  use  mean  squares,  normalized  correlation, 
and  mutual  information  to  demonstrate  the  accuracy  and  efficiency  of  the  multiscale  registration 
technique. 


1  Introduction 

Often  in  image  processing,  images  must,  be  spatially  aligned  in  order  to  perform  quantitative  analyses  of 
the  images.  The  process  of  aligning  images  taken,  for  example,  at  different  times,  from  different  imaging 
devices,  or  from  different  perspectives,  is  called  image  registration.  More  precisely,  image  registration  is 
the  process  of  determining  the  optimal  spatial  transform  that  maps  one  image  to  another.  Typically,  two 
images  are  taken  as  input,  and  the  registration  process  is  then  the  optimization  problem  which  determines 
the  geometric  mapping  that  brings  one  image  into  spatial  alignment  with  the  other  image.  In  practice,  the 
particular  type  of  transformation  as  well  as  the  notion  of  optimal  will  depend  on  the  specific  application. 

Examples  of  applications  in  which  image  registration  is  particularly  important  include  astro-  and  geo¬ 
physics,  computer  vision,  remote  sensing,  and  medicine.  In  this  paper,  we  will  focus  on  medical  image 
registration.  Image  registration  plays  an  important  role  in  the  analysis  of  medical  images.  For  example, 
images  taken  from  different  sensors  often  contain  complementary  information.  By  bringing  the  two 
images  into  alignment  so  that  anatomical  features  of  one  modality  can  be  detected  in  the  other  modality, 
the  information  from  the  different  modalities  can  be  combined.  In  neurosurgery,  for  example,  tumors  are 
typically  identified  and  diagnosed  using  magnetic  resonance  images  (MRI) ,  but  stereotaxy  technology  (the 
use  of  surgical  instruments  to  reach  specified  points)  generally  uses  computed  tomography  (CT)  images. 
Registration  of  these  modalities  allows  the  transfer  of  coordinates  of  tumors  from  the  MRI  images  to 
the  CT  images.  See  [13]  for  a  discussion  of  the  applications  of  multi-modality  imaging  to  problems  in 
neurosurgery.  As  another  example,  medical  image  data  acquired  prior  to  diagnosis  can  be  compared  with 
data  acquired  during  or  after  treatment  to  determine  the  effectiveness  of  the  treatment.  To  compare 
images  taken  at  different  times,  however,  the  images  must  first  be  brought  into  spatial  alignment  with 
one  another  so  that  actual  differences  in  the  data  can  be  distinguished  from  differences  that  result  from 
the  image  acquisition  process. 

In  the  context  of  medical  imaging,  the  goal  of  the  registration  process  is  to  remove  artificial  differences 
in  the  images  introduced  by  patient  movement,  differences  in  imaging  devices,  etc.,  but  at  the  same 
time,  to  retain  real  differences  due  to  actual  variations  of  the  objects.  Medical  images,  however,  often 
contain  significant  levels  of  noise  due  to  instrumentation  imperfections,  data  acquisition  techniques,  image 
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reconstruction  methods,  transmission  and/or  compression  errors,  and  other  factors.  Although  numerous 
successful  image  registration  techniques  have  been  published,  we  will  see  that  ordinary  image  registration 
algorithms  can  fail  to  produce  meaningful  results  when  one  or  both  of  the  images  to  be  registered  contains 
significant  levels  of  noise. 

Since  noise  is  generally  present  in  digital  images,  image  denoising  is  a  fundamental  problem  in  image 
processing.  Indeed,  numerous  approaches  to  image  denoising  have  been  presented.  Thus  a  simple  solution 
to  the  problem  of  image  registration  in  the  presence  of  noise  would  be  to  first  apply  a  denoising  algorithm 
to  the  noisy  image(s),  and  then  use  existing  image  registration  techniques  to  register  the  denoised  images. 
However,  common  denoising  algorithms,  most  notably  spatial  filtering  algorithms,  have  the  disadvantage 
that  while  they  are  successful  in  removing  noise,  they  often  remove  edges  as  well.  Additionally,  most 
denoising  procedures  require  a  priori  knowledge  of  the  noise  level,  variance,  and/or  model,  while  this 
information  is  typically  not  known  in  practice.  For  these  and  other  reasons,  we  will  demonstrate  that 
ordinary  image  registration  of  noisy  images  fails  to  produce  acceptable  results  even  when  classical  denois¬ 
ing  algorithms  are  applied  to  the  noisy  images  prior  to  registration  (for  significantly  high  levels  of  noise). 
Thus  we  seek  a  technique  that  enables  successful  image  registration  when  one  or  both  of  the  images  to 
be  registered  contains  noise. 

In  practice,  we  can  consider  an  image  /  of  consisting  of  coarse  and  fine  scales.  The  general  shape  and 
main  features  of  an  image  are  considered  the  coarse  scales,  and  details  and  textures,  such  as  noise,  are  the 
fine  scales  of  the  image.  Separating  the  coarse  and  fine  scales  of  an  image,  therefore,  is  an  effective  tool 
in  denoising.  Indeed,  several  denoising  algorithms  have  been  proposed  using  separation  of  the  coarse  and 
fine  scales  of  an  image,  most  notably  [18],  [17],  [10],  and  [19].  The  method  presented  in  [19]  presents  a 
multiscale  technique  in  which  an  image  /  is  decomposed  in  a  hierarchical  expansion  /  ~  £ jUj,  where  the 
iifs  (called  the  components  of  /  relative  to  the  decomposition)  resolve  edges  of  /  with  increasing  scales. 
More  precisely,  for  small  k,  the  sum  Ej'uj  is  a  coarse  representation  of  the  image  /,  and  as  k  increases, 
the  sum  captures  more  and  more  detail  (and  hence,  noise)  of  the  image. 

In  this  paper,  we  present  a  multiscale  image  registration  technique  based  on  the  multiscale  decom¬ 
position  of  [19]  that  is  particularly  effective  when  one  or  both  of  the  images  to  be  registered  contains 
significant  levels  of  noise.  Since  the  hierarchical  expansion  /  ~  £ jUj  decomposes  the  image  /  into  com¬ 
ponents  which  contain  increasingly  fine  scales,  we  expect  a  component-wise  registration  algorithm  to 
produce  accurate  results  for  noisy  images.  That  is,  given  a  noisy  image  /,  for  small  values  of  k,  the  com¬ 
ponent  E jUj  retains  the  general  shape  of  the  image  /  while  removing  the  details  and  noise  of  the  image. 
Thus  if  we  wish  to  register  the  image  /  with  another  image,  say  g  we  expect  that  registration  of  the 
components  Ej-Uj  with  g  will  provide  an  accurate  estimation  of  the  actual  transformation  that  brings  the 
two  images  into  spatial  alignment  with  one  another,  for  sufficiently  small  values  of  k.  Similarly,  if  both  / 
and  g  are  noisy,  we  expect  decomposing  both  images  and  performing  component- wise  registrations  should 
accurately  estimate  the  optimal  transformation.  We  will  demonstrate  that  multiscale  image  registration 
enables  successful  image  registration  for  images  that  contain  levels  of  noise  significantly  higher  than  the 
levels  at  which  ordinary  registration  fails. 

This  paper  is  organized  in  the  following  way.  In  Section  3,  we  discuss  the  image  registration  problem 
and  review  standard  image  registration  techniques.  In  Section  4,  we  present  the  problem  of  image 
registration  in  the  presence  of  noise,  and  illustrate  the  failure  of  current  techniques  when  one  or  both  of  the 
images  to  be  registered  contains  high  levels  of  noise.  We  also  briefly  discuss  classical  denoising  techniques, 
and  illustrate  the  failure  of  ordinary  image  registration  of  noisy  images  even  when  the  images  are  denoised 
prior  to  registration.  In  Section  6,  we  review  the  multiscale  image  decomposition  of  [19],  and  illustrate 
the  results  of  the  hierarchical  multiscale  decomposition  obtained  upon  applying  the  algorithm  to  noisy 
images.  In  Section  7,  we  present  image  registration  techniques  based  upon  the  multiscale  decomposition 
described  in  Section  6,  and  in  Section  8,  we  present  our  registration  results.  Finally,  in  Section  9,  we 
discuss  computational  aspects  of  our  registration  algorithm. 
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3  The  registration  problem 

Given  a  fixed  and  moving  image,  the  registration  problem  is  the  process  of  finding  an  optimal  transfor¬ 
mation  that  brings  the  moving  image  into  spatial  alignment  with  the  fixed  image.  While  this  problem  is 
easy  to  state,  it  is  difficult  to  solve.  The  main  source  of  difficulty  is  that  the  problem  is  ill-posed,  which 
means,  for  example,  that  the  problem  may  not  have  a  unique  solution.  Additionally,  the  notion  of  opti¬ 
mality  may  vary  for  each  application:  for  example,  some  applications  may  require  consideration  only  of 
rigid  transformations,  while  other  applications  require  non-rigid  transformations.  Finally,  computation 
time  and  data  storage  constraints  place  limitations  on  the  complexity  of  models  that  can  be  used  for 
describing  the  problem. 

3.1  The  mathematical  setting 

A  two-dimensional  gray-scale  image  /  is  a  mapping  which  assigns  to  every  point  x  E  C  R2  a  gray 
value  f{x)  (called  the  intensity  value  of  the  image  at  the  point  x).  We  will  consider  images  as  elements 
of  the  space  L2(R2).  Color  images  can  be  defined,  for  example,  in  terms  of  vector-valued  functions 
f  =  (/i :  / 2 :  /a)  representing  the  RGB-color  scales.  For  the  medical  imaging  applications  that  we  are 
interested  in,  images  are  in  fact  given  in  terms  of  discrete  data,  and  the  function  /  must  be  obtained  via 
interpolation.  We  will  not  discuss  this  construction  here,  but,  assume  that  an  interpolation  method  has 
been  chosen  and  fixed. 

Image  registration  is  typically  necessary  when  two  images  are  essentially  of  the  same  object,  but  the 
images  are  not  spatially  aligned.  This  occurs,  for  example,  when  the  images  are  taken  at  different  times, 
from  different  perspectives,  or  from  different  imaging  devices.  The  basic  input  data  to  the  registration 
process  are  two  images:  one  is  defined  as  the  fixed  image  f(x)  and  the  other  as  the  moving  image  m(x). 
The  goal  is  then  to  find  a  transformation  0  such  that  the  fixed  image  /(.?:)  and  the  deformed  moving 
image  m$(x)  :=  m(4>(x))  are  similar.  To  solve  this  problem  in  a  mathematical  way,  the  term  similar 
needs  to  be  defined  in  an  appropriate  fashion.  For  example,  if  the  images  to  be  registered  are  taken 
from  different  devices,  there  may  not  be  a  correspondence  between  the  intensities  f(x)  and  m^x)  for 
an  optimal  <j>.  Additionally,  we  may  consider  measures  of  similarity  between  the  images  which  are  not 
related  to  the  intensities.  Thus  the  registration  problem  necessarily  involves  a  discussion  of  the  distance 
measures,  or  metrics,  used  to  compare  images.  The  general  problem  can  then  be  stated  as  follows: 

Given  a  distance  measure  D  :  (L2(R2))2  — >  R  and  two  images  f(x),m(x)  E  L2(R2),  the  solution  of  the 
registration  problem  is  given  by  the  following  minimization  problem: 

4>  —  argmin  D (/,  m (3.1) 

^:R2^R2 

In  many  applications,  the  set  of  allowable  transformations  to  be  considered  in  the  minimization 
problem  (3.1)  is  restricted  to  a  strict  subset  of  the  set  of  all  maps  if  :  R2  — >  R2.  For  example,  we  may 
require  the  transformation  cb  to  be  smooth,  or  we  may  impose  specific  parametric  requirements,  such  as 
requiring  the  transformation  to  be  rigid,  affine,  polynomial,  or  spline. 

3.2  Landmark-based  registration 

Landmark-based  registration  is  an  image  registration  technique  which  is  based  on  a  finite  set  of  image 
features.  The  problem  is  then  to  determine  the  transformation  such  that  for  a  finite  set  of  features,  any 
feature  of  the  moving  image  is  mapped  onto  the  corresponding  features  of  the  fixed  image.  More  precisely, 
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let  F (f,j)  and  F j  —  1, . . .  ,m  be  given  features  of  the  fixed  and  moving  images,  respectively.  The 
solution  4>  of  the  registration  problem  is  then  a  map  cfc  :  R2  — *  JR2  such  that 

F(/,j)  =  4>{E{m,j)),  j  =  (3.2) 

For  a  more  general  notion  of  landmark-based  registration,  we  define  the  following  distance  measure: 

m 

Dlm(0)  :=  £  ||F(/,j)  -  mm-,  j)\\l  (3-3) 

i=i 

where  ||  •  ||;  denotes  a  norm  on  the  landmark,  or  feature,  space.  For  example,  if  the  features  are  locations 
of  points,  then  ||  •  ||;  =  ||  •  |[R2.  We  can  then  restate  (3.2)  as  the  minimization  problem  in  which  the 
solution  cf> :  R2  — *  K2  of  the  registration  problem  is  given  by: 

(j>  =  argmin  D LM(V-')-  (3-4) 

V>:R2-»R2 

To  solve  this  minimization  problem,  the  transformation  is  chosen  to  either  be  an  element  of  an  n- 
dimensional  space  spanned,  for  example,  by  polynomials,  splines,  or  wavelets,  or  it  is  required  to  be 
smooth  in  some  sense.  In  the  first  case,  the  features  to  be  mapped  are  the  locations  of  a  number  of 
user-supplied  landmarks.  Let  Xk,  k  =  1, . . .  n  be  the  basis  functions  of  the  space.  Then  the  minimization 
of 

m  m 

DLM(^)  :=  J]  ||F (f,j)  -  =  DLM(«/>)  :=  £  ||F (f,j)  -  ^(F(m,  j)||R2 

j= i  j=i 

can  be  obtained  upon  expanding  <f>  =  {4>i ,  fa)  in  terms  of  the  basis  functions  and  solving  the  resulting 
least  squares  problems. 

In  the  case  in  which  we  require  the  transformation  <j>  to  be  smooth,  we  introduce  a  functional  which 
imposes  smoothness  restrictions  on  the  transformation.  That  is,  we  look  for  a  transformation  (f>  which 
interpolates  the  features  F (f,  j)  and  F(m,j),  and  which  is  smooth  in  some  sense.  Such  a  transformation 
is  called  a  minimal  norm  solution,  and  it  turns  out  (see  [8])  that  the  solution  can  be  expressed  in  terms 
of  radial  basis  functions. 

Landmark-based  registration  is  simple  to  implement  and  the  numerical  solution  requires  only  the 
solution  of  a  linear  system  of  equations.  However,  the  main  drawback  of  the  landmark-based  approach 
is  that  the  registration  process  depends  on  the  location  of  the  landmarks.  As  the  detection  and  math¬ 
ematical  characterization  of  landmarks  (for  example,  anatomical  landmarks  in  medical  images)  is  not 
fully  automated,  the  landmarks  must  be  user-supplied,  and  this  can  be  a  time-consuming  and  difficult 
process,  even  for  a  medical  expert;  see,  for  example,  [16].  Additionally,  landmark-based  registration  does 
not  always  results  in  a  physically  meaningful  registration.  See  [11],  pp.  44,  for  a  simple  example  of  a 
situation  in  which  landmark-based  registration  fails  to  produce  meaningful  results. 


3.3  Principal  axes-based  registration 

Principal-axes  image  registration  is  based  on  the  idea  of  landmark-based  registration,  but  it  uses  features 
that  can  be  automatically  detected.  The.se  features  are  constructed  as  follows.  For  an  image  /  :  JR2  — +  M, 
and  a  function  g  :  R2  — >  M,  we  define  the  expectation  value  of  g  with  respect  to  /  by 


%(<?)  == 


Jr2  g(x)f{x)  dx 
fX2  g{x)  dx 


If  u  :  R2  — y  MmX71,  we  set  E/(u)  :=  (E /[%,&]  €  Rmxn. 


(3.5) 
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The  center  of  an  image  /  is  defined  by 


C  f.=  £/[:>:]  €  R2 

and  the  covariance  by 


(3.6) 


Co vj  :=  Ej  [(x  -  Cj)(x  -  C/)T]  G  R2x2  (3.7) 

Given  fixed  and  moving  images  /(./;)  and  m(x),  the  centers  cj  and  c.m  and  eigendecompositions  of 
the  covariance  matrices  Cov/  and  Covm  are  used  as  the  features  F,,  and  the  registration  problem  is  to 
compute  <j) :  R2  — »  R2  such  that  Fj(rn(^>))  =  F ,  (/)  for  the  features  F,;. 

This  method  is  described  in  detail  in  [lj.  The  prinicipal-axes  method  of  image  registration  has  the 
advantages  that  it  is  computationally  fast  and  simple  and  requires  few  registration  parameters,  but  has 
the  disadvantages  that  it  is  not  suitable  for  images  of  multiple  modalities  and  that  the  solutions  may  be 
ambiguous.  In  particular,  the  prinicipal-axes  based  method  cannot  distinguish  between  images  with  the 
same  center  and  covariance,  even  though  images  with  very  different  structure  and  orientation  may  have 
the  same  center  and/or  covariance. 

3.4  Optimal  parametric  registration 

In  this  section,  we  present  methods  of  registration  which  are  based  on  the  minimization  (or  maximization) 
of  some  distance  measure,  or  metric,  D.  The  transformation  <j>  is  restricted  to  some  parameterized  space, 
and  the  registration  can  be  ovtained  by  minimizing  (or  maximizing)  the  distance  D  over  the  space.  In 
particular,  we  will  discuss  metrics  based  on  intensity,  correlation,  and  mutual  information.  Given  a  metric 
D,  a  fixed  image  /,  and  a  moving  image  m,  optimal  parametric  registration  is  the  problem  of  finding 
a  transformation  <j)  in  some  pre-specified  parameterizable  space  such  that  D(/,  m(<f>))  is  minimized  (or 
maximized  in  certain  cases).  Examples  of  commonly  used  parameterizable  spaces  in  image  registration  are 
polynomial  and  spline  spaces.  We  will  primarily  be  interested  in  rigid  and  affine  linear  transformations. 
An  affine  linear  map  is  a  map  of  the  form  cj>(x)  =  Aa;  +  b,  A  €  R2x2,  det  A  >  0,  b  G  R2.  Such  a  map  allows 
rotations,  translations,  scales,  and  shears  of  the  coordinates.  A  translation  (or  rigid)  transformation  is  a 
special  case  of  an  affine  transformation  which  allows  only  rotations  and  translations  of  the  coordinates, 
and  in  this  case,  the  matrix  A  is  required  to  be  orthogonal  with  determinant  1.  Optimal  parametric 
registration  is  probably  the  most  commonly  used  image  registration  technique. 

3.4.1  Choosing  an  optimizer 

To  minimize  D  (/,  m(0)),  we  must  choose  an  optimization  technique.  That  is,  an  optimal  parametric 
registration  technique  is  described  by  a  metric  to  be  minimized  (or  maximized)  and  an  optimizer  which 
controls  the  minimization  (or  maximization).  The  implementation  of  the  registration  algorithm  works 
in  the  following  way.  At  each  iteration,  the  distance  D  between  the  two  images  is  computed.  An  affine 
transformation  is  then  applied  to  the  moving  image,  and  the  distance  between  the  images  is  recomputed. 
This  process  continues  until  the  distance  is  minimized  (or  maximized) . 

The  optimizer  controls  this  algorithm.  At  each  stage,  the  optimizer  determines  the  parameters  of 
the  transformation  that  will  be  applied  to  the  moving  image.  Examples  of  commonly  used  optimizers 
include  gradient  descent  and  regular  step  gradient  descent.  Gradient  descent  optimization  advances  the 
parameters  of  the  transformation  in  the  direction  of  the  gradient,  where  the  step  size  is  governed  by  a 
user-specified  learning  rate.  Regular  step  gradient  descent  optimization  advances  the  parameters  of  the 
transformation  in  the  direction  of  the  gradient  where  a  bipartition  scheme  is  used  to  compute  the  step 
size. 
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3.4.2  Mean  squares  metric 

The  mean  squares  metric  computes  the  mean-squared  pixel-wise  difference  in  intensity  between  two 
images  f  and  m: 

1  N 

MS(/,  m)  :=  —  -  mi )2>  (3-8) 

N  i= 1 

where  N  is  the  total  number  of  pixels  considered,  /,•  is  the  *th  pixel  of  image  /,  and  rrq  is  the  *th 
pixel  of  image  m.  Note  that  the  optimum  value  of  the  mean  squares  metric  is  0,  and  poor  matches 
between  the  images  /  and  m  results  in  large  values  of  MS (f,m).  This  metric  has  the  advantage  that  it 
is  computationally  simple,  but,  it  is  based  on  the  assumption  that  pixels  in  one  image  should  have  the 
same  intensity  as  (spatially)  corresponding  pixels  in  the  second  image.  Thus,  the  mean  squares  metric  is 
restricted  to  images  of  the  same  modality. 

3.4.3  Normalized  correlation  metric 

The  normalized  correlation  metric  computes  pixel-wise  cross-correlation  and  normalizes  it  by  the  square 
root  of  the  autocorrelation  function: 


E  (fi-mi) 

NC (/,  rn)  :=  - 1  x  ■  ,i=1  ■  ,  (3.9) 

/  N  N 

y  i=i  2=i 

where  N,  /*,  and  m,  are  as  defined  for  the  mean  squares  metric.  The  -1  factor  in  (3.9)  causes  the  optimum 
value  of  the  metric  to  occur  when  the  minimum  is  reached.  Thus  the  optimal  value  of  the  normalized 
correlation  metric  is  -1.  As  with  the  mean  squares  metric,  the  normalized  correlation  metric  is  restricted 
to  images  of  the  same  modality. 

3.4.4  Viola- Wells  mutual  information  metric 

Mutual  information  is  an  information-theoretic  approach  to  image  registration  that  was  proposed  inde¬ 
pendently  by  Viola  and  Wells  [21]  and  Collignon  et  al  [4]  in  1995.  The  idea  is  that  mutual  information 
computes  the  amount  of  information  that  one  random  variable  (here,  image  intensity)  gives  about  another 
random  variable  (here,  intensity  values  of  another  image).  More  precisely,  given  a  fixed  image  f(x)  and  a 
moving  image  m(x),  we  wish  to  compute  the  transformation  6  which  maximizes  the  mutual  information 

(p  =  arg  max  I(f(x),m(<f>(x))).  (3.10) 

The  maximization  of  mutual  information  criterion  assumes  that  the  statistical  dependence  between  cor¬ 
responding  image  intensity  values  is  maximized  when  the  images  are  geometrically  aligned. 

The  mutual  information  I (f(x),m(cf>(x))  is  defined  in  terms  of  entropy,  where  we  consider  i  as  a 
random  variable  over  coordinate  locations  in  the  coordinate  system  of  the  fixed  image.  Let  //,(■)  denote 
the  entropy  of  a  random  variable:  h(x)  :=  —  f  p(x)lnp(x)  dx.  where  p(x)  is  the  probability  density 
function  of  the  random  variable  x.  Note  that  it  is  not  clear  how  to  construct  p(x);  we  will  discuss 
methods  for  estimating  the  probability  densities.  The  joint  entropy  of  two  random  variables  x  and  y  is 
given  by  h(x,  y)  =  —  J  p(x,  y)  lnp(x,  y)dxdx,  where  p(x,  y )  is  the  joint  probability  density  function  of  the 
random  variables  x  and  y.  Entropy  can  be  considered  as  a  measure  of  the  uncertainty  or  complexity  of 
a  random  variable. 
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If  x  and  y  are  independent,  then  p(x,y)  =  p(x)p(y),  so  h(x,y )  =  h(x )  +  h{y).  However,  if  there  is 
any  dependency  (as  would  be  the  case  if  x  and  y  are  intensity  values  of  images  of  the  same  object),  then 
h{x,y)  <  h(x)  +  h(y).  The  difference  is  defined  to  be  mutual  information: 

I(f(x),  m{4>{x)))  =  h(f(x))  +  h(m{<j>(x))  -  h{f{x),  (3.11) 

The  terms  in  (3.11)  can  be  interpreted  in  the  following  way.  The  first  term  h(f(x))  is  the  entropy  of 
the  fixed  image  and  is  independent  of  the  transformation  <j>.  The  second  term  h(m(<j>(: :;;)))  is  the  entropy 
of  m(<j>(x )),  so  maximization  of  mutual  information  encourages  transformations  <j>  for  which  m.(<f>(x)  has  a 
high  level  of  complexity  or  uncertainty.  The  third  term  -h{f{x),m(<j>(x)))  is  the  negative  joint  entropy  of 
f(x)  and  m(<j)(x)),  so  maximization  of  mutual  information  is  related  to  minimization  of  the  joint  entropy 
of  f(x)  and  m(<j)(x)).  Taken  together,  the  second  and  third  terms  identify  transformations  <j>  that  find 
complexity  in  the  images  and  explain  it  well.  In  [15],  the  authors  present  a  detailed  overview  of  mutual 
information  based  registration. 

Mutual  information  has  the  following  properties.  Let  u(x)  and  v(x)  denote  two  images. 

1.  I (u(x),v(x))  =  l(v(x).  u(x  j)  Mutual  informations  is  symmetric.  Although  this  is  true  theoretically, 
in  practice,  it  is  not  always  the  case  that  we  obtain  the  same  transformation  upon  registering  u(x) 
with  v(x)  and  v(x)  with  u(x).  This  is  a  consequence  of  numerical  implementation  methods. 

2.  I(u(x),  u(x))  =  h(u(x)).  The  information  an  image  contains  about  itself  is  equal  to  the  entropy  of 
the  image. 

3.  I(ii(x),  u(x))  <  h(u(x))  I(u(x),v(x))  <  h(v(x))  The  information  that  the  images  contain  about,  each 
other  can  not  be  greater  than  the  information  contained  in  the  individual  images. 

4.  I(u(x),  v(x))  >  0 

5.  I(u(.x),  v(.t))  =  0  if  and  only  if  u(x)  and  v(x)  are  independent.  That  is,  if  the  images  u(x)  and  v(x) 
are  independent,  no  information  about  one  image  is  gained  when  the  other  image  is  known. 


The  entropies  in  Equation  (3.11)  are  defined  in  terms  of  integrals  over  the  probability  densities  as¬ 
sociated  with  the  images  f{x)  and  m(x).  However,  in  a  typical  medical  image  registration  problem,  the 
probability  densities  are  not  directly  accessible,  and  thus  must  be  estimated  from  the  image  data.  Parzen 
windowing,  described  in  [5],  is  a  commonly  used  technique  for  estimating  the  densities,  and  it  is  the 
one  used  by  Viola  and  Wells  in  [21].  In  this  method,  continuous  density  functions  are  constructed  by  a 
super-position  of  kernel  functions  K(  )  centered  the  elements  of  a  sample  of  intensities  taken  from  the 
image.  That  is,  the  estimation  of  the  probability  density  p(z)  is  given  by: 


p(x)  =  P*(z) 


(3.12) 


where  Ns  is  the  number  of  spatial  samples  in  S  and  K  is  an  appropriately  chosen  kernel  function.  The 
kernel  function  K  must  be  smooth,  symmetric,  have  zero  mean,  and  integrate  to  1.  Examples  of  suitable 
candidates  for  K  include  the  Gaussian  density  and  the  Cauchy  density.  In  [21],  Viola  and  Wells  use  a 
Gaussian  density  function  with  standard  deviation  a  to  estimate  the  probability  density  functions.  The 
optimal  value  of  a  depends  on  the  particular  images  to  be  registered. 

Upon  estimating  the  probability  densities  using  the  Parzen  windowing  technique,  the  entropy  integral 
h(z)  =  —  / p(z)  ln(p(z))  dz  must  be  evaluated.  This  integral  is  difficult,  or  impossible,  to  evaluate 
analytically,  so  it  must  be  approximated  as  a  sample  mean: 
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h(z)  =  ~~  E  ln(p*(^))>  (3-13) 

R  ZjSR 

where  R  is  a  second  sample  of  intensities  taken  from  the  image.  That  is,  two  separate  intensity  samples 
S  and  R  are  taken  from  the  image.  The  first  is  used  to  estimate  the  probability  density,  and  the  second 
is  used  to  approximate  the  entropy. 

The  main  advantage  of  the  mutual  information  measure  is  that  it  is  generally  applicable  for  multi¬ 
modality  registration,  whereas  intensity-based  measures  are  typically  not  applicable  for  multi-modality 
registration.  Mutual  information  registration  has  been  used  successfully  for  a  number  of  difficult,  applica¬ 
tions.  Most  notable,  mutual  information  has  been  shown  to  be  highly  accurate  for  MRI-CT  registration. 
See,  for  example,  [9],  [14],  and  [20],  for  a  discussion  of  these  applications. 

3.5  Non-parametric  image  registration 

All  of  the  image  registration  techniques  that  we  have  discussed  so  far  have  been  based  on  certain  parame¬ 
ters.  For  example,  either  the  transformation  <j>  can  be  expanded  in  terms  of  basis  functions  that  span  a 
specified  finite-dimensional  space,  or  the  registration  is  controlled  by  a  specified  set  of  external  features. 
Non-parametric  techniques  do  not  restrict  the  transformation  to  a  parameterizable  set.  Given  two  images, 
a  fixed  image  f(x )  and  a  moving  image  m(x),  we  seek  a  transformation  0  such  that  m(cj>(x))  is  similar  to 
fix')  in  a  certain  sense.  Upon  defining  a  suitable  distance  measure  D,  the  registration  problem  is  then  to 
minimize  the  distance  between  and  f{x).  However,  a  direct  minimization  is  often  not  possible 

in  the  non-parametric  case.  For  example,  the  problem  is  ill-posed:  small  changes  in  the  input  data  may 
lead  to  large  changes  in  the  output.  Additionally,  the  solution  is  not  unique.  Given  these  constraints,  a 
stable  numerical  implementation  is  often  not  possible.  To  circumvent  these  problems,  a  regularizing,  or 
smoothing,  term  S  is  introduced,  and  the  registration  problem  becomes  the  minimization  of  the  distance 
between  m(<f>(x))  and  f(x)  plus  a  smoothing  term  That  is,  the  registration  is  based  on  a  regularized 
minimization  of  the  distance  between  the  images. 

In  the  discussion  of  non-parametric  image  registration,  the  transformation  :  R2  — >  K2  is  split  into 
the  trivial  identity  part  and  the  deformation  or  displacement  part  u,  i.e. 

4>{x)  =  x  -  u(x).  (3.14) 

Upon  decomposing  <f>  in  this  way,  we  have  m(4>(x))  =  m(x  —  u(x))  :=  mu(x).  Given  a  distance  D  and 
a  smoother  S,  the  elastic  registration  problem  is  then  the  minimization  of  D(  f(x),  mu(x))+aS(u),  where 
a  €  M  is  a  positive  regularizing  parameter. 

The  choice  of  smoother  S  typically  depends  on  the  particular  application.  Examples  of  non-parametric 
image  registration  techniques  include  elastic  registration  [3],  fluid  registration  [2],  and  diffusion  registra¬ 
tion  [7].  Elastic  registration  uses  linear  elasticity  theory  to  model  the  deformation  of  an  elastic  body. 
In  this  case,  the  regularizing  term  S (u)  is  the  linearized  elastic  potential  of  the  displacement  u.  In  fluid 
registration,  the  regularization  is  based  on  the  linearized  elastic  potential  of  the  time  derivative  of  u. 
Finally,  diffusion  registration  uses  regularization  based  on  spatial  derivatives  of  the  displacement. 

3.6  Remarks 

In  this  section,  we  presented  a  brief  overview  of  the  major  image  registration  techniques  currently  used 
in  medical  image  registration.  In  practice,  the  best  registration  method  for  a  given  set  of  images  will 
depend  on  the  particular  features  of  the  images  themselves.  However,  numerous  comparison  studies  which 
compare  the  accuracy  and  performance  of  different  image  registration  techniques  for  various  applications 
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have  been  presented.  The  most  extensive  of  these  is  [22],  which  originally  consisted  of  a  comparison  of 
16  different  methods,  but  has  since  been  substantially  expanded. 


4  Registration  in  the  presence  of  noise 

In  this  section,  we  study  the  effect  of  noise  on  image  registration,  and  we  determine  the  approximate 
noise  level  at  which  registration  fails.  Consider  the  brain  proton  density  slice  images  shown  in  the  figure 
below.  The  image  on  the  right  is  the  result  of  translating  the  image  on  the  left  by  13  ram  in  X  and  17 
mm  in  Y.  Let  I  denote  the  original  image,  and  let  T  denote  the  translated  image. 


Image  I 

Brain  Proton  Density  Slice 


Translated  Image  T 


Figure  1:  Original  image  I  and  translated  image  T 

Initially,  we  will  consider  the  registration  problem  in  which  one  of  the  images  (here,  the  fixed  image) 
is  noisy.  We  will  add  increasing  levels  of  noise  to  the  image  I  and  register  the  the  non-noisy  translated 
image  T  with  the  noisy  images.  Our  goal  is  to  determine  the  approximate  noise  levels  at  which  various 
image  registration  techniques  fail,  and  to  develop  an  algorithm  that  will  enable  registration  beyond  these 
levels.  Eventually,  we  will  also  apply  our  techniques  to  the  case  in  which  both  the  fixed  and  moving 
images  contain  significant  levels  of  noise.  Before  we  present  these  results,  we  discuss  the  notion  of  noise 
in  more  mathematical  detail. 

4.1  Noise 

Digital  images  are  often  degraded  by  random  noise.  In  imaging,  the  term  noise  refers  to  random  fluc¬ 
tuations  in  intensity  values  that  occur  during  image  capture,  transmission,  or  processing,  and  that  may 
distort  the  information  given  by  the  image.  Image  noise  is  not  part  of  the  ideal  signal  and  may  be 
caused  by  a  wide  range  of  sources,  such  as  detector  sensitivity,  environmental  radiation,  transmission 
errors,  discretization  effects,  etc.  Noise  is  generally  classified  as  either  independent  noise,  or  noise  which 
is  dependent  on  the  image  data. 
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Independent  noise  can  often  be  described  by  an  additive  noise  model,  in  which  the  observed  image  / 
is  the  sum  of  the  true  image  s  and  the  noise  n: 

f(x)  -  s(x) +  n(x).  (4.1) 

Within  this  framework  of  additive  noise,  the  noise  n{x)  is  commonly  modeled  by  Gaussian  noise  of  mean 
m  and  variance  v. 

A  multiplicative  noise  model  describes  noise  that  is  dependent  on  the  image  data.  This  is  often 
referred  to  as  speckle  noise. 


f(x)  =  s(x)  +  s(a:)n(x)  =  s(z)(l  +  n(  x))  (4.2) 

In  this  case,  n(x )  is  uniformly  distributed  random  noise  with  mean  m  and  variance  v. 

Impulse  noise,  or  salt  and  pepper  noise,  is  noise  that  resembles  salt  and  pepper  granules  randomly 
distributed  over  the  image.  Impulse  noise  is  typically  defined  by  the  following  model.  Again  we  let  s(x) 
denote  the  actual  image,  and  f(x)  denote  the  observed  image. 


/(*)  = 


\ti{x) 


with  probability  1  —  S 
with  probability  <5, 


(4.3) 


where  r](x)  is  an  identically  distributed,  independent  random  process.  With  this  model,  an  arbitrary 
pixel  x  €  fl  c  M2  is  affected  by  noise  with  probability  5,  and  not  affected  with  probability  1  —  <5.  We 
will  refer  to  5  as  the  impulse  noise  density,  as  adding  impulse  noise  of  density  S  to  an  image  f(x)  affects 
approximately  5  ■  siz e(/)  pixels.  The  random  process  rj(x)  is  typically  such  that  the  corrupted  pixels  are 
either  set  to  the  maximum  value,  have  single  bits  flipped  over,  or  are  set  alternatively  to  zero  or  to  the 
maximum  value.  This  last  case  results  in  a  "‘salt  and  pepper"’  appearance.  Note  that  unaffected  pixels 
always  remain  unchanged. 

In  Figure  2,  we  add  additive  Gaussian  noise  of  mean  0  and  variance  0.2,  multiplicative  speckle  noise 
of  mean  0  and  variance  0.2,  and  impulse  noise  of  density  0.2  to  the  brain  proton  density  slice  image  I. 

In  this  paper,  we  will  study  the  problem  of  image  registration  in  the  presence  of  high  levels  of  impulse 
noise,  although  we  believe  that  our  solution  is  applicable  to  registration  in  the  presence  of  other  forms 
of  noise,  as  well.  We  add  impulse  noise  of  increasing  densities  5  to  the  brain  proton  density  slice  image 
I,  and  register  the  (non-noisy)  translated  image  T  with  the  noisy  images.  Let  Ig  denote  the  image  I  with 
added  impulse  noise  of  density  6.  In  3,  we  illustrate  the  noisy  images  for  increasing  values  of  6. 


4.2  Registration  results 

For  each  5  in  3,  we  register  T  with  lg  using  various  registration  methods.  Recall  that  the  image  T  is 
the  result  of  translating  the  original  image  I  13  units  in  X  and  17  units  in  Y,  and  that  la  is  the  result 
of  adding  uniform  impulse  noise  of  density  S  to  the  image  I.  Since  T  is  a  rigid  transformation  of  I, 
we  will  restrict  the  registration  process  to  linear  transformations,  i.e.  we  will  consider  optimal  linear 
registrations.  The  optimal  transformation  <j>  produced  by  the  optimal  linear  registration  process  will 
consist  of  two  parameters,  namely  X-  and  Y-translation  values.  We  will  let  rpx  and  rpy  denote  the  X- 
and  Y-translation  parameters,  respectively,  of  the  optimal  transformation  <fi.  For  comparison  purposes, 
we  will  perform  the  optimal  linear  registration  using  the  mean  squares  metric,  normalized  correlation 
metric,  and  mutual  information  metric. 

We  use  the  following  parameters  for  the  registration  algorithms.  For  the  mean  squares  and  normalized 
correlation  registration  algorithms,  we  use  the  regular  step  gradient  descent  optimizer.  Due  to  the 
stochastic  nature  of  the  metric  computation  in  the  mutual  information  algorithm,  the  regular  step  gradient 
descent  optimizer  does  not  work  well  in  the  case  of  mutual  information.  Instead,  we  use  the  gradient 
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Additive  Gaussian 
noise  of  mean  0 
and  variance  0.2 


Multiplicative  speckle 
noise  of  mean  0 
and  variance  0.2 


Impulse  noise  of 
density  0.2 


Figure  2:  In  this  figure,  we  illustrate  the  addition  of  various  types  of  noise  to  the  image  I. 


Original  Image 


8=0.10 


8=0.20 


8=0.30 


Figure  3:  In  this  figure,  we  illustrate  the  addition  of  impulse  noise  of  increasing  densities  5 
to  the  image  I. 
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Mean  squares 

Normalized  correlation 

Mutual  information 

6 

4>x 

4>y 

n 

4>x 

<t> Y 

n 

<fat 

4>y 

n 

0 

12.99 

17.00 

18 

13.01 

17.00 

18 

12.75 

17.03 

200 

0.10 

12.99 

17.01 

28 

12.99 

17.01 

20 

12.83 

16.88 

200 

0.20 

13.03 

16.98 

17 

13.04 

16.98 

19 

12.98 

16.64 

200 

0.30 

12.97 

17.03 

28 

13.02 

17.02 

11 

13.02 

17.02 

200 

0.40 

18.89 

7.16 

15 

8.05 

1.30 

13 

11.08 

9.72 

200 

0.50 

2.16 

7.06 

19 

9.09 

2.18 

8 

9.72 

7.12 

200 

0.60 

29.81 

3.19 

40 

4.08 

0.24 

7 

4.57 

5.17 

200 

0.70 

2.08 

1.14 

13 

3.11 

2.13 

12 

3.08 

2.86 

200 

Table  1:  In  this  table,  we  present  the  results  obtained  upon  registering  the  translated  image  T  with  the 
noisy  image  1^ ,  where  8  is  the  impulse  noise  density.  We  let  4>x  and  4>y  denote  the  X-  and  Y-translation 
values  of  the  optimal  transformation  <f>  produced  by  the  registration  algorithm,  and  we  denote  by  n  the 
number  of  iterations  until  convergence.  Recall  that  the  actual  translation  values  are  13  units  in  X  and 
17  units  in  Y. 


descent  optimizer  with  a  user-specified  learning  rate  of  20.0.  Finally,  we  set  the  maximum  number  of 
iterations  for  each  algorithm  to  200.  As  we  will  see,  mean  squares  and  normalized  correlation  registrations 
typically  converge  very  quickly  to  the  optimum  value.  Mutual  information,  on  the  other  hand,  often  does 
not  ever  actually  reach  the  true  optimal  solution,  but  instead  oscillates  within  one  or  two  pixels  of  the 
optimal  solution  (generally  after  100-150  iterations).  By  reducing  the  learning  rate,  we  can  increase  the 
likelihood  of  convergence,  but  this  increases  the  computation  time  significantly  without  improving  the 
accuracy  of  the  solution. 

For  each  of  these  three  registration  algorithms,  and  for  each  5  we  record  the  X-  and  Y-translation 
parameters,  denoted  (f>\  and  <j> y,  respectively,  of  the  optimal  transformation  0  produced  by  the  registration 
process.  We  also  record  the  number  of  iterations  n  until  convergence.  The  results  are  illustrated  in  Table 
1.  Recall  that  the  actual  translation  values  are  13  units  in  X  and  17  units  in  Y.  We  also  record  the 
number  of  iterations  until  convergence,  which  we  denote  by  n. 

The  results  presented  in  Table  1  indicate  that  optimal  linear  registration  in  the  presence  of  impulse 
noise  fails  when  the  impulse  noise  density  in  the  fixed  image  reaches  approximately  0.40. 

5  Denoising 

5.1  Denoising  techniques 

In  this  section,  we  discuss  various  denoising  techniques.  Image  denoising  is  a  fundamental  problem 
in  image  processing,  and  there  has  been  much  research  and  progress  on  the  subject.  As  our  primary 
interest  in  this  paper  is  the  problem  of  image  registration  of  noisy  images,  and  not  denoising,  we  do  not 
focus  on  the  general  problem  of  image  denoising,  but  instead  present  a  few  of  the  most  common  and 
computationally  simple  denoising  techniques.  We  will  then  apply  these  techniques  to  one  of  our  noisy 
images,  and  study  the  effect  of  denoising  on  the  image  registration  techniques.  In  particular,  in  Section 
4,  we  saw  that  ordinary  optimal  linear  registration  of  noisy  images  failed  when  the  impulse  noise  density 
was  greater  than  0.40.  In  this  section,  we  determine  whether  or  not  denoising  prior  to  registration  enables 
successful  registration  of  noisy  images  for  which  registration  failed  previously. 

Spatial  filtering  is  the  traditional  approach  to  removing  noise  from  images.  Spatial  filters  use  the 
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assumption  that  noise  occupies  the  higher  regions  of  the  frequency  spectrum,  and  thus  attenuate  high 
spatial  frequencies.  Filtering  is  a  neighborhood  process,  in  which  the  value  of  a  given  pixel  in  the 
filtered  image  is  computed  by  applying  some  algorithm  to  the  pixel  values  in  a  neighborhood  of  the  given 
pixel.  Typical  implementations  of  spatial  filters  include  mean  filtering,  median  filtering,  and  Gaussian 
smoothing.  Mean  filtering  computes  the  value  of  each  output  pixel  by  computing  the  statistical  mean  of 
the  neighborhood  of  the  corresponding  input  pixel.  Thus,  applying  a  mean  filter  to  a  noisy  image  reduces 
the  amount  of  variation  in  gray-level  intensity  between  pixels.  Although  this  filter  is  computationally 
easy  to  implement,  it  is  sensitive  to  the  presence  of  outliers.  Median  filtering,  which  computes  the  value 
of  each  output  pixel  by  computing  the  statistical  median  of  the  neighborhood  of  the  corresponding  input 
pixel,  is  more  robust  to  the  presence  of  outliers,  and  is  thus  commonly  used  for  removing  impulse  noise 
from  images.  Convolution  with  a  Gaussian  kernel  is  another  commonly  used  spatial  filtering  technique. 

In  Figure  4,  we  illustrate  the  effect  of  applying  a  mean,  median,  and  Gaussian  convolution  filter  to 
the  noisy  image  I0.70,  the  brain  proton  density  slice  image  with  impulse  noise  of  density  0.70. 


Noisy  image 


Mean  filter 


Median  filter  Gaussian  convolution 


Figure  4:  In  this  figure,  we  illustrate  the  results  of  applying  three  different  denoising  filters  to  the  brain 
proton  density  slice  image  with  impulse  noise  of  density  0.70. 
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Mean  squares 

Normalized  correlation 

Mutual  information 

Denoising  Technique 

</>x 

<j)y 

n 

</>x 

<t>  Y 

n 

< fix 

4>y 

n 

Mean  filtering 

31.83 

1.15 

46 

16.88 

1.11 

29 

5.39 

5.30 

200 

Median  filtering 

18.87 

1.26 

31 

2.38 

6.90 

34 

4.39 

4.06 

200 

Gaussian  convolution 

18.86 

-0.76 

31 

2.19 

0.25 

11 

7.38 

7.37 

200 

Table  2:  In  this  table,  we  present  the  results  obtained  upon  registering  the  translated  image  T  with  the 
denoised  images  obtained  upon  applying  median,  mean,  and  Gaussian  convolution  filters  to  the  noisy 
image  I0.70.  We  let  <f>\  and  4>v  denote  the  X-  and  Y-translat.ion  values  of  the  optimal  transformation  <j> 
produced  by  the  registration  algorithm,  and  we  denote  by  n  the  number  of  iterations  until  convergence. 
Recall  that  the  actual  translation  values  are  13  units  in  X  and  17  units  in  Y. 

5.2  Registration  results  after  denoising 

In  this  section,  we  register  the  translated  image  T  with  the  denoised  images  illustrated  in  Figure  4.  As 
in  Section  4,  we  use  mean  squares,  normalized  correlation,  and  mutual  information  optimal  linear  reg¬ 
istration.  For  each  registration  method,  we  let  <f>  denote  the  optimal  transformation  produced  by  the 
registration  algorithm,  and  we  let  4>x  and  <f>y  the  X-  and  Y-translation  parameters  of  the  optimal  trans¬ 
formation  cf>.  We  denote  by  n  the  number  of  iterations  of  each  registration  algorithm  until  convergence. 
We  record  the  results  in  Table  2.  The  moving  image  in  each  case  is  the  translated  image  T;  recall  that 
the  actual  translation  values  are  13  in  X  and  17  in  Y. 

The  results  presented  in  Table  2  indicate  that  the  application  of  common  denoising  techniques  prior 
to  registration  does  not  enable  successful  registration  of  the  noisy  image  I0.70  with  the  translated  image 
T.  While  the  denoising  algorithms  presented  here  are  reasonably  successful  in  removing  noise,  they 
also  remove  a  significant  amount  of  detail  from  the  images  and  blur  edges.  Since  denoising  prior  to 
registration  fails  to  give  successful  results,  we  present  an  image  registration  technique  based  on  a  multiscale 
decomposition  of  the  image(s)  to  be  registered. 


6  Multiscale  decomposition 

6.1  The  hierarchical  decomposition 

In  this  section,  we  present  the  multiscale  image  representation  using  hierarchical  (BV,  L2)  decompositions 
of  [19].  Consider  an  image  /  €  L2(0).  Define  the  J-functional  J (/,  A)  as  follows: 

J (/,•*)  :=  |nf  {A|M|22  +  ||w||Bv}, 

U+V=f 

where  A  >  0  is  a  scaling  parameter  that  separates  the  L2  and  BV  terms.  This  functional  J(/,  A)  was 
introduced  in  the  context  of  image  processing  by  Rudin,  Osher,  and  Fatemi  [18].  They  suggest  the 
following.  Let  \u\,  ua]  denote  the  minimizer  of  J(/,  A).  The  BV  component,  u\  captures  the  coarse 
features  of  the  image  /,  while  the  L2  component,  v\  captures  the  finer  features  of  /  such  as  noise.  This 
model  is  effective  in  denoising  images  while  preserving  edges,  though  it  requires  prior  knowledge  on  the 
noise  scaling  A. 

Tadmor,  Nezzar,  and  Vese  propose  an  alternative  point  of  view  in  which  the  minimization  of  J(/,  A)  is 
interpreted  as  a  decomposition  /  =  u\  +  v\,  where  u\  extracts  the  edges  of  /  and  v\  extracts  the  textures 
of  /.  This  interpretation  depends  on  the  scale  A,  since  "‘texture"’  at  scale  A  consists  of  edges  when  viewed 
under  a  refined  scale  (2A,  for  example).  Then  we  decompose  v\  as  follows: 

v\  =  u2\  +  v2\,  where  [i^a^a]  =  arginf  J('<;a,  2A). 

U-\~V=V  \ 
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Thus  we  obtain  a  two-scale  representation  of  /  given  by  /  =  u\+u 2\.  Continuing  this  process,  we  obtain 
a  hierarchical  multiscale  decomposition  of  /,  as  follows.  Starting  with  an  initial  scale  A  =  Ao,  we  obtain 
an  initial  decomposition  of  the  image  /: 


f  =  u0  +  v0,  [u0,  «o]  =  arginf  J(/;  A0). 

U+V=f 

We  then  refine  this  decomposition  to  obtain 

Vj  =  uj+ 1  +  vj+i,  [uj+i,  Vj+i]  =  arginf  3(vj,  A02?+1),  j  =  0,1,.. .. 

U-\-V~Vj 


After  k  steps  of  this  process,  we  have  the  following  hierarchical  decomposition  of  /: 
/  =  u0  +  v0 

=  Mo  +  Ul  +  Wl 

=  U0+Ui  +  U2+V2 


—  U0  +  Ul  +  ■  ■  ■  Uk  +  Vk 

Thus  we  obtain  a  multiscale  image  decomposition  /  ~  uo  +  U\  +  . . .  +  Uk,  with  a  residual  Vfc.  As  k 
increases,  the  u k  components  resolve  edges  with  increasing  scales  A*  =  \)2k ■ 

6.2  Implementation 

6.2.1  Initialization 

As  described  in  ||,  the  initial  scale  Ao  should  capture  the  smallest  oscillatory  scale  in  /,  given  by 

2W  ^  ll/llw-'.oo  < 

However,  in  practice,  we  may  not  be  able  to  determine  the  size  of  ||/||w-i,~,  so  we  determine  the  initial 
choice  of  Ao  experimentally.  For  the  applications  presented  in  this  paper,  we  will  use  Ao  =  0.01  and 
A  j  =  A02T 


6.2.2  Numerical  discretization 

We  follow  the  numerical  algorithm  of  Tadmor,  Nezzar,  and  Vese  for  the  construction  of  our  hierarchical 
decomposition.  In  each  step,  we  use  finite-difference  discretization  of  the  Euler-Lagrange  equations  as¬ 
sociated  with  the  Aj+i)  to  obtain  the  next,  term,  Uj+ 1,  in  the  decomposition  of  the  image  /.  The 
Euler-Lagrange  equation  associated  with  the  minimization  of  J (/,  A)  is 

-  5xdiv(iwxr)  =  /. 

with  Neumann  boundary  conditions. 

k 

We  thus  obtain  an  expansion  /  ~  E  where  the  Uj  are  constructed  as  approximate  solutions  of  the 

3= 0 

recursive  relation  given  by  the  following  elliptic  PDE: 

U3+ 1  -  2A7WdiV(F^tl)  =  -2AJdiV(^7T)- 

Note  that  J(/,  A)  contains  a  singularity  when  |Vwa|  =  0.  To  remove  this  singularity,  we  replace  J(/,  A) 
by  the  regularized  functional 

Jc(/,  A)  :=  inf  {AIM&  +  fn  ^ t2  +  |Vm|2  dx  dy}, 

U+V=J 

and  at  each  step,  we  find  the  minimizer  ii\  of  Je.  The  Euler-Lagrange  equation  for  the  regularized  J( 
functional  is 
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u\ 


Adiv( 


Viix 


v/^+ Tv^IF 


)  =  /  in  0, 


with  Neumann  boundary  conditions. 

To  numerically  implement  the  method,  we  cover  the  domain  0  with  a  grid  (a;*  :=  th,  :=  jh),  and 
discretize  the  elliptic  PDE  given  above  using  the  forward,  backward,  and  centered  divided  differences. 
Let  D+,  D_,  and  Do  denote  the  forward,  backward,  and  centered  divided  differences,  respectively.  The 
discretized  PDE  is  as  follows: 


fi,j  +  2A^ 
>-» 

=  f*,j  +  2S2 


+  ^D_„[  i 

2A  Vl  v/£2  +  (D0l«i.j)!*  +  (D+t(«i,iP 


V'c2  +  (D+x-u,.JV  +  (Do,/«i.ipD+;cUi’:'i 

D +yUij\ 


ui  +  l.j  ui,j 


Uj—  l ,  j 


+ 


y/c?+(r>+*Ui,j)2+(tooyUi,j)*  \/  c2+{D_3,«i1j)2+(Doy 

1  f  _ i j  +  1  ~Ui>3 _  _  Ui,3~ ui,j~l  _ 1 

2*IV£2+(Do*«i,i)2+(D+*«ij)2  v/£2  +  (Do-“‘.i-»)2  +  (D-»«i.i)2 


To  solve  the  discrete  regularized  Euler-Lagrange  equations,  we  use  the  Gauss-Siedel  iterative  method  to 
obtain: 


,,«+ 1 


—  +  2ft2 


J  ^  2ft2  1  \/£2+(D+,«L)2+(Dos«"j.)2  v/e2+(D_a;UV.)2  +  (Dos,<_I,J.)2j 

. n  _ n+1  n+1  n 


+ 


2ft2  1  v/e2  +  (D0x<J)2+(D+v<,,.)2  <ScM(D0xuZj_1)2+(D-v<,j')2 


To  satisfy  the  Neumann  boundary  conditions,  we  first  reflect  /  outside  fi  by  adding  grid  lines  on  all 
sides  of  fl.  As  the  initial  condition,  we  set  uJP  =  We  iterate  this  numerical  scheme  for  n  =  0, 1, . . .  N 
until  ||un=“  —  —  1  [  |  is  less  than  some  preassigned  value  so  that  u”“  is  an  accurate  approximation  of 

the  fixed  point  steady  solution  u\. 

Finally,  we  denote  the  final  solution  u\  :=  }tJ.  To  obtain  the  hierarchical  multiscale  decompo¬ 
sition,  we  reiterate  this  process,  each  time  updating  /  and  A  in  the  following  way: 


fnew  <  t current  U\, 

^new  <  2A current- 


That  is,  at  each  step,  we  apply  the  J(/ current  —  u\,  2 A)  minimization  to  the  residual  f current  ~  uiambda 
of  the  previous  step.  Taking  A  j  =  Ao2J,  we  obtain  after  k  steps  a  hierarchical  multiscale  decomposition 
f  =  U\0  +  Uiambda x  +  •  • .  +  u\k  +  v\k ,  where  we  write  uxj  —  Uj.  We  call  the  uj,  j  =  1,  2, . . . ,  k  the 
components  of  /  and  the  v *  the  residuals. 

6.3  Decomposition  of  a  noisy  image 

Consider  the  image  I0.70  shown  below.  This  is  the  brain  proton  density  slice  image  I  with  impulse  noise 
of  density  0.70. 

We  apply  the  hierarchical  multiscale  decomposition  to  this  noisy  image,  using  the  following  parame¬ 
ters: 

•  k  =  12  hierarchical  steps 

•  A0  =  0.01 
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Figure  5:  Noisy  Image  Iq  70 


•  A  j  =  Ao2^ 

•  e  =  0.001 

•  n  =  10 

•  h.=  1 

The  figures  below  illustrate  the  components  maj  and  the  residuals  v^j  for  this  decomposition.  Note  that 
in  each  hierarchical  step,  an  additional  amount  of  texture  is  seen  in  the  components.  Further,  the  noise 
is  not  seen  in  the  first  few  components,  while  most  of  the  texture  is  kept,  and  the  noise  only  reappears 
as  the  refined  scales  reach  the  same  scales  as  the  noise  itself.  Eventually,  we  will  use  this  multiscale 
decomposition  to  register  the  noisy  image  I0.70  with  the  translated  image  T. 

7  Multiscale  Registration 

Consider  the  brain  proton  density  slice  I  with  added  impulse  noise  of  density  S.  Recall  that  registration  of 
the  the  translated  image  T  with  the  noisy  image  failed  when  5  >  0.40  using  the  mean  squares,  normalized 
correlation,  and  mutual  information  registration  methods.  Moreover,  registration  using  these  classical 
methods  failed  even  after  denoising  the  noisy  image  using  various  standard  denoising  techniques.  In 
this  section,  we  present  new  methods  for  image  registration  that  enable  successful  registration  of  the 
translated  image  T  with  the  noisy  images  0  for  values  of  the  noise  density  <5  significantly  greater  than  the 
levels  at  which  classical  registration  and  registration  after  denoising  fails.  These  registration  techniques 
will  be  based  on  the  hierarchical  multiscale  decomposition  described  in  Section  6. 

Consider  two  images  A  and  B,  and  suppose  that  we  want  to  register  image  B  with  image  A.  Suppose 
that  one  or  both  of  the  images  contains  a  significant  amount  of  noise.  If  only  one  of  the  images  is  noisy, 
call  the  noisy  image  A.  We  propose  the  following  multiscale  registration  method.  First,  we  apply  the 
multiscale  hierarchical  decomposition  to  both  images.  Let  m  denote  the  number  of  hierarchical  steps 

k 

used  for  the  multiscale  decompositions.  For  ease  of  notation,  given  an  image  /,  we  let  C *,(/)  :=  'uAfc 

i- 0 

denote  the  fcth  component  of  the  image  f,k  —  0, 1, ... ,  m  —  1  obtained  as  in  Section  6.  Thus  Cfc(A)  will 
denote  the  kth  component  of  the  image  A,  and  C*,(B)  will  denote  the  kth  component  of  image  B. 

We  will  present  two  separate  algorithms;  in  the  first,  we  consider  registration  of  image  B  with  compo¬ 
nents  of  image  A,  and  in  the  second,  we  consider  registration  of  components  of  image  B  with  components 
of  image  A. 
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Figure  6:  Multiscale  decomposition  of  the  noisy  image  I0.70 
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7.1 

In  our  first  multiscale  registration  algorithm,  we  register  image  B  with  the  kth  component  of  A,  for 
k  =  0, 1, ... ,  m  —  1.  This  is  illustrated  by  the  following  schematic: 

Let  fa  denote  the  optimal  transformation  produced  by  the  registration  algorithm  upon  registering  B 
with  Ca-(A),  k  =  0, 1, . . . ,  m— 1.  Recall  that  Co(A)  contains  only  the  coarsest  scales  of  the  image  A,  and  as 
k  increases,  (A  (A)  contains  increasing  levels  of  detail  (and  hence,  noise)  of  the  image  A.  Thus  we  expect 
that  registration  of  image  B  with  Cfc(A)  should  give  an  improvement  compared  to  ordinary  registration 
for  the  first  few  values  of  k.  As  k  increases,  however,  we  expect  that  eventually  the  component  Cfc(A) 
will  become  too  noisy  to  give  successful  registration. 

Upon  determining  the  transformations  fa  via  a  suitable  registration  algorithm,  we  have  several  options 
for  defining  the  optimal  transformation  $  that  should  bring  the  image  B  into  spatial  alignment  with  the 
image  A.  The  first  option  would  be  to  define  <f>  as  an  ordinary  average  of  the  fa: 

~  7/1—1 

*  ==  -  £  fa-  (7.1) 

,U  k= 0 

However,  as  we  mentioned  above,  we  do  not  expect  the  transformation  fa  to  be  an  accurate  estimation 
of  the  actual  optimal  transformation  $  for  large  values  of  k.  Thus,  a  second  option  would  be  to  define  $ 
as  a  weighted  average  of  the  fa. 


^  rn— 1 

$  :=  —  yZ  0>k<t>ki  (7.2) 

m  f— ' 
k= 0 

where  the  weights  a *  are  appropriately  chosen  non-negative  real  numbers.  We  will  consider  these  defini¬ 
tions  in  more  detail  when  we  present  our  results  in  Section  8. 


7.2 

In  our  second  multiscale  registration  algorithm,  we  register  the  ktb  component  of  image  B  with  the  kth 
component  of  image  A,  for  k  =  0, 1,  2, ...  m  —  1,  as  illustrated  in  the  following  schematic: 

Let  fa  denote  the  optimal  transformation  produced  by  the  registration  algorithm  upon  registering 
Cfc(B)  with  Cfc(A),  k  =  0, 1, ... , m  —  1.  As  before,  we  expect  that  registration  of  Cfc(B)  with  C^(A)  should 
give  an  improvement  compared  to  ordinary  registration  for  the  first  few  values  of  k.  As  k  increases, 
however,  we  expect  that  eventually  the  components  Cfc(A)  and  C/;(B)  will  become  too  noisy  to  give 
successful  registration.  Since  this  algorithm  considers  components  of  both  images,  we  expect  that  it  will 
be  particularly  successful  in  the  case  in  which  both  images  are  noisy. 

Finally,  we  define  the  optimal  transformation  T  that  should  bring  image  B  into  spatial  alignment 
with  image  A  using  either  an  ordinary  average: 


or  a  weighted  average: 


m—  1 


*  :=  -  £  fa; 

m  f— ' 


k=0 


^  m—  1 

^  :=  —  Y]  htpk, 

rn  a — y 


k~0 


where  the  weights  bk  are  appropriately  chosen  non-negative  real  numbers. 


(7.3) 


(7.4) 
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Mutual  information 

Mean  squares 

Normalized  correlation 

Fixed  Image 

4>x 

<t>Y 

n 

<t>x 

<t>Y 

n 

<f>x 

4>y 

n 

Noisy  image  I0.70 

4.57 

5.17 

200 

2.08 

1.14 

7 

4.08 

0.24 

7 

O 

b- 

d 

HH 

O 

O 

12.65 

16.36 

200 

3.08 

1.11 

12 

3.11 

0.17 

9 

Cl(Io.7o) 

12.69 

16.78 

200 

2.08 

3.08 

14 

2.13 

2.12 

12 

C2(Io.7o) 

12.56 

16.79 

200 

2.11 

3.08 

14 

2.14 

3.11 

15 

C3(Io.7o) 

12.53 

16.76 

200 

3.08 

2.11 

14 

3.11 

2.14 

7 

O 

t— 1 
O 

O 

12.48 

16.76 

200 

24.88 

1.16 

36 

18.86 

1.18 

30 

C5  (Io.7o) 

12.46 

16.78 

200 

40.80 

1.07 

52 

0.21 

1.18 

11 

Ce(Io.7o) 

12.43 

16.80 

200 

28.86 

0.15 

46 

27.84 

2.19 

42 

C7(Io.7o) 

12.43 

16.79 

200 

-2.87 

4.11 

15 

0.18 

3.14 

12 

C8(Io.7o) 

12.43 

16.74 

200 

25.89 

3.12 

40 

-1.84 

4.12 

14 

C9(Io.7o) 

9.33 

9.41 

200 

6.05 

4.12 

12 

7.99 

2.08 

16 

Clo(Io.7o) 

8.44 

8.32 

200 

-3.92 

8.12 

21 

4.09 

3.15 

16 

Cll(Io.7o) 

6.96 

6.46 

200 

8.97 

6.13 

13 

3.65 

1.17 

27 

Table  3:  This  table  illustrates  the  registration  results  upon  registering  the  translated  image  T  with  the 
kth  component  C*.-  (Io.ro)  of  the  noisy  image  I0.70  obtained  via  the  multiscale  decomposition  discussed 
in  Section  6.  Here,  we  use  m  —  12  hierarchical  steps  to  decompose  the  noisy  image,  so  we  perform 
m  =  12  registration  simulations.  The  transformation  parameters  4>x  and  cpy  are  the  X-  and  Y-translation 
parameters  of  the  optimal  transformation  <p  produced  by  the  registration  algorithm.  Recall  that  the  actual 
translation  values  are  13  in  X  and  17  in  Y.  The  moving  image  in  all  simulations  is  the  translated  image 
T.  The  moving  image  in  each  case  is  the  translated  image  T. 

8  Multiscale  registration  results 

8.1  Noisy  fixed  image 

In  this  section,  we  use  the  multiscale  registration  technique  described  in  Section  7  to  register  the  translated 
(non-noisy)  image  T  with  the  noisy  image  I0.70-  Recall  that  I0.70  is  the  image  obtained  upon  adding 
impulse  noise  of  density  0.70  to  the  brain  proton  density  slice  image  I.  As  before,  let  Ct(Io.7o)  denote  the 
kth  component  in  the  multiscale  decomposition  of  Io.7o,  for  k  =  0, 1, . .  .  m,  obtained  as  in  Section  7.  We 
perform  the  multiscale  decomposition  using  m  —  12  hierarchical  steps,  Ao  =  0.01,  and  Xj  —  Xq2L  In  the 
table  below,  we  present  the  results  of  m  =  12  registration  simulations,  obtained  upon  registering  T  with 
Cfc(Io.7o),  k  =  0, 1, ... ,  11.  For  each  registration,  we  let  4>  denote  the  optimal  transformation  produced 
by  the  registration  algorithm,  and  we  let  cj>x  and  <j>y  the  X-  and  Y-translation  parameters  of  the  optimal 
transformation  <j>.  The  moving  image  in  each  registration  is  the  translated  image  T.  For  reference,  we 
also  include  in  the  first  line  of  Table  3  the  parameters  obtained  using  ordinary  registration. 

To  provide  an  estimation  of  the  optimal  transformation  $  that  should  bring  the  translated  image  T 
into  alignment  with  the  noisy  image  I0.70,  we  first  compute  an  ordinary  average  as  described  in  Equation 
(7.1).  Letting  4>x  and  'Ey  denote  the  corresponding  X-  and  Y-translation  parameters  of  the  average  4>, 
we  obtain: 

It  is  clear  from  the  results  presented  in  Table  3,  however,  that  a  weighted  average  as  constructed 
in  Equation  7.2  is  more  appropriate.  In  particular,  we  see  that  the  parameters  produced  by  mutual 
information  registration  are  clustered  around  12.5  units  in  X  and  16.8  units  in  Y  for  k  =  0, 1, . . .  8,  but 
then  are  significantly  different  for  the  remaining  values  of  k.  We  expected  that  the  multiscale  registration 
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Mutual  information 

Mean  squares 

Normalized  correlation 

11.45 

11.6 

6.29 

$Y 

14.02 

3.11 

2.05 

Table  4:  This  table  gives  the  translation  parameters  4>x  and  obtained  by  computing  the  ordinary 
average  of  the  multiscale  registration  transformations  <fit  for  k  —  0, 1, ...  11,  as  described  in  Equation 
(7.1).  Recall  that  the  actual  translation  values  are  13  in  X  and  13  in  Y. 


Mutual  information 

12.52 

4?y 

16.73 

Table  5:  In  this  table,  we  compute  the  weighted  average  of  the  multiscale  registration  transformations 
<f>k  for  k  —  0,1,...  11  for  mutual  information  registration.  We  thus  obtain  the  corresponding  X-  and 
Y-translation  parameters  of  the  transformation  4>  defined  as  in  Equation  (7.1).  Recall  that  the  actual 
translation  values  are  13  in  X  and  17  in  Y. 


results  would  be  an  accurate  approximation  of  the  actual  transformation  fI>  for  small  values  of  k,  but 
then  would  deviate  as  k  became  sufficiently  large,  because  as  k  becomes  large,  increasing  scales  of  detail 
(and  hence,  noise)  appear  in  the  component  Ck  ■  Thus,  even  without  knowing  the  actual  values  of  the  X- 
and  Y-translations,  it  makes  sense  to  consider  in  the  weighted  average  for  mutual  information  only  the 
parameters  corresponding  to  the  first  9  registrations  ( k  =  0, 1, ...  8).  We  define  the  weights  a*.  =  1  for 
k  =  0, 1, . .  .8  and  =  0  for  k  >  9.  In  Table  5,  we  give  the  resulting  X-  and  Y-translation  parameters. 
Since  the  actual  values  are  13  in  X  and  17  in  Y,  we  see  that  multiscale  mutual  information  registration 
produced  very  accurate  results  in  this  case,  and  indeed  is  a  significant  improvement  compared  to  ordinary 
registration  as  well  as  to  classical  denoising  followed  by  registration. 

Note  that  we  could  use  more  sophisticated  techniques  for  determining  the  weights  a For  example, 
we  could  perform  a  statistical  analysis  on  the  parameters  given  by  the  <j>k,  and  compute  the  weights  a*, 
according  to  the  mean  and  standard  deviation.  Finally,  we.  see  that  the  parameters  produced  by  the 
mean  squares  and  normalized  correlation  registrations,  however,  do  not  cluster  around  a  single  value.  In 
these  cases,  the  registration  algorithm  did  not  produce  a  meaningful  result. 

Next,  we  provide  the  results  obtained  upon  registering  the  multiscale  components  of  the  translated 
image  T  with  the  multiscale  components  of  the  noisy  image  Io.ro-  Let  Cfc(T)  and  C/t(Io.7o)  denote  the 
multiscale  components  of  T  and  I0.70,  respectively,  obtained  via  the  multiscale  decomposition  presented 
in  Section  6.  As  before,  we  use  m  =  12  hierarchical  steps,  Ao  =  0.01,  and  A j  =  XoV  to  perform  the 
decomposition.  In  Table  6,  we  present  the  results  of  m  =  12  registration  simulations,  obtained  upon 
registering  Cfc(T)  with  C*(Io.7o),  k  —  0,1,...,  11.  For  each  registration,  we  let  ip  denote  the  optimal 
transformation  produced  by  the  registration  algorithm,  and  we  let  ipx  and  ij!Y  denote  the  X-  and  Y- 
translation  parameters  of  the  optimal  transformation  ip.  For  reference,  we  also  include  in  the  first  line  of 
Table  6  the  parameters  obtained  using  ordinary  registration. 

The  results  in  Table  6  indicate  that  a  weighted  average  as  constructed  in  Equation  (7.4)  is  the  most 
appropriate  method  for  estimating  the  optimal  transform  parameters  and  Recall  that  we  expect 
the  multiscale  registration  algorithm  to  work  well  for  the  first  few  values  of  k.  and  then  to  eventually 
become  less  accurate  as  increasing  levels  of  detail  and  noise  appear  in  the  kth  component  C*,.  Indeed, 
for  mutual  information  registration,  we  see  that  the  parameters  ipx  and  ipy  are  clustered  together  for 
k  =  0, 1, ...  8,  and  then  deviate  for  k  >  9,  so  we  set  the  weights  =  1  for  k  =  0, 1, ...  8  and  a*,  =  0  for 
k  >  9.  For  both  mean  squares  and  normalized  correlation,  we  set  the  weights  bk  =  1  for  k  —  0, 1  and 
bk  =  0  for  k  >  1.  In  Table  7,  we  present  the  X-  and  Y-translation  values  corresponding  to  the  weighted 
average  '3>. 
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Mutual  information 

Mean  squares 

Normalized  correlation 

Fixed 

Image 

Moving 

Image 

ipx 

1pY 

n 

i>x 

Ipy 

n 

ipx 

Ipy 

n 

1-0.70 

T 

4.57 

5.18 

200 

2.08 

1.14 

13 

4.08 

0.24 

7 

Co(Io.7o) 

Co(T) 

12.69 

16.66 

200 

12.29 

17.72 

21 

12.96 

17.08 

17 

Cl(Io.7o) 

Ci(T) 

12.67 

16.87 

200 

13.70 

17.75 

38 

12.99 

17.67 

16 

C2(Io.7o) 

C2(T) 

12.59 

16.86 

200 

20.77 

5.20 

38 

16.84 

4.31 

17 

Cs(Io.7o) 

Cs(T) 

12.55 

16.82 

200 

3.19 

0.31 

12 

4.20 

4.23 

11 

C4(Io.7o) 

C4(T) 

12.52 

16.83 

200 

2.20 

2.24 

12 

26.74 

5.18 

36 

(l-0.7o) 

C5(T) 

12.51 

16.84 

200 

31.65 

2.23 

42 

14.90 

6.27 

58 

Ce(Io.7o) 

Ce(T) 

12.49 

16.87 

200 

30.69 

6.16 

50 

19.87 

4.29 

28 

C7(Io.7o) 

c7(T) 

12.48 

16.85 

200 

33.64 

3.16 

46 

29.64 

3.32 

37 

o 

00 

»— 1 
© 

“'I 

O 

C8(T) 

12.53 

16.71 

200 

28.81 

3.22 

45 

1.26 

1.29 

11 

Cg(Io.7o) 

C9(T) 

9.26 

9.36 

200 

2.13 

3.13 

11 

17.93 

3.21 

71 

Clo(Io.7o) 

Cio(T) 

8.80 

8.61 

200 

2.12 

3.12 

11 

32.63 

3.14 

39 

Cll(Io.7o) 

Cn(T) 

6.95 

6.34 

200 

34.74 

2.10 

42 

4.13 

5.08 

17 

Table  6:  This  table  illustrates  the  registration  results  upon  registering  the  kth  multiscale  component 
Cfc(T)  of  the  translated  image  T  with  the  kth  multiscale  component  Cfc(Io.7o)  of  the  noisy  image  I0.70 
obtained  via  the  multiscale  decomposition  discussed  in  Section  6.  Here,  we  use  m  =  12  hierarchical 
steps  to  decompose  the  noisy  image,  so  we  perform  m  =  12  registration  simulations.  The  transformation 
parameters  ipx  and  ipY  are  the  X-  and  Y-translation  parameters  of  the  optimal  transformation  ip  produced 
by  the  registration  algorithm.  Recall  that  the  actual  translation  values  are  13  in  X  and  17  in  Y. 


Mutual  information 

Mean  squares 

Normalized  correlation 

S'x 

12.56 

12.99 

12.98 

4/y 

16.82 

17.74 

17.37 

Table  7:  This  table  gives  the  translation  parameters  d'x  and  d'y  obtained  by  computing  the  weighted 
average  of  the  multiscale  registration  transformations  </;fc  for  k  =  0,1, ...11,  as  described  in  Equation 
(7.4).  Recall  that  the  actual  translation  values  are  13  in  X  and  17  in  Y. 
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Since  the  actual  translation  values  are  13  in  X  and  17  in  Y,  we  see  that  the  component- wise  multiscale 
registration  of  the  translated  image  T  with  the  noisy  image  I0.70  produces  very  accurate  results  for  each 
of  the  three  methods  presented  here  (mean  squares,  normalized  correlation,  and  mutual  information). 

8.2  Noisy  fixed  and  moving  images 

In  this  section,  we  consider  the  registration  problem  in  which  both  the  fixed  and  moving  images  are  noisy. 
Consider  the  noisy  images  I0.40  and  T0.40,  where  T,  as  before,  is  the  result  of  translating  I  13  units  in  X 
and  17  units  in  Y,  and  As  denotes  the  image  obtained  by  adding  impulse  noise  of  density  8  to  the  image 
A. 


Figure  8:  Original  image  I  and  translated  image  T  with  impulse  noise  of  density  5  =  0.40 

Before  applying  our  multiscale  registration  algorithm,  we  attempt  to  register  the  noisy  translated 
image  T0.40  with  the  noisy  image  I0.40  using  the  three  registration  methods  mean  squares,  normalized 
correlation,  and  mutual  information.  For  each  registration  method,  we  denote  by  <f>  the  optimal  transfor¬ 
mation  produced  by  the  registration,  and  we  denote  by  <f>x  and  4>y  the  X-  and  Y-translation  parameters 
of  the  optimal  transformation  <j>.  As  is  seen  in  the  table  below,  registration  of  the  noisy  images  fails. 

Since  ordinary  registration  of  the  noisy  images  fails,  we  register  the  images  using  our  multiscale 
registration  technique.  First,  we  perform  the  multiscale  decomposition  discussed  in  Section  6  to  both  of 
the  noisy  images,  again  using  k  =  12  hierarchical  steps,  initial  scale  Ao  =  0.01,  and  A3  =  ‘P  ■  An-  Let 
Cfc(Io.4o)  and  Cfc(To.4o)  denote  the  kth  component  in  the  multiscale  decomposition  of  I0.40  and  T0.40, 
respectively.  Since  both  of  the  images  are  noisy,  we  register  the  fcth  component  Cfc(T0,4o)  with  the 
kth  component  Cfc(Io.4o)-  For  each  registration  simulation,  we  denote  by  ■</;  the  optimal  transformation 
produced  by  the  registration  algorithm,  and  we  denote  by  ijjx  and  ipy  the  corresponding  X-  and  Y- 
translation  parameters  of  the  optimal  transformation  ?/>. 

The  multiscale  registration  results  in  Table  9  indicate  that  a  weighted  average  as  constructed  in 
Equation  (7.2)  is  the  most  appropriate  method  for  estimating  the  optimal  transform  parameters  fhx  and 
Ty-  For  mutual  information,  we  set  the  weights  6*,  =  1  for  k  =  0, 1, . . .  6  and  6*.  =  0  for  k  >7.  For  mean 
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Registration  Method 

<Px. 

<Py 

Mean  Squares 

11.02 

7.04 

Normalized  Correlation 

3.05 

0.99 

Mutual  Information 

5.03 

2.54 

Table  8:  This  table  illustrates  the  results  upon  registering  the  noisy  translated  image  T0.40  with  the  noisy 
image  I0.40  using  ordinary  mean  squares,  normalized  correlation,  and  mutual  information  registration 
techniques.  For  each  registration,  we  let  <p  denote  the  optimal  transformation  produced  by  the  registration 
algorithm,  and  we  let  <px  and  <}>y  denote  the  corresponding  X-  and  Y-translation  parameters  of  the  optimal 
transformation  <p.  Recall  that  the  actual  translation  values  are  13  in  X  and  17  in  Y. 


Mutual  information 

Mean  squares 

Normalized  correlation 

Fixed 

Image 

Moving 

Image 

V’x 

i’y 

ipx 

1pY 

ipx 

1pY 

I0.40 

To. 40 

5.03 

2.54 

11.02 

7.04 

3.05 

0.99 

Co(Io.4o) 

Co(To.4o) 

13.06 

16.92 

13.05 

16.92 

13.05 

16.92 

Cl(I0.40) 

Ci(To.40) 

13.05 

16.93 

13.02 

16.22 

13.06 

16.92 

C2(Io.4o) 

C2(To.4o) 

13.03 

16.93 

8.11 

5.29 

13.02 

16.27 

C3  (I0.40) 

C3(T0.4o) 

13.02 

16.94 

5.40 

12.19 

13.02 

16.25 

O 

O 

O 

C4(T0,4o) 

13.02 

16.94 

2.20 

8.00 

2.23 

5.09 

C5(Io.4o) 

C5(T0.4o) 

13.01 

16.93 

26.76 

1.21 

1.17 

7.00 

Ce(Io.4o) 

C6(T0.4o) 

12.99 

16.81 

23.83 

4.11 

1.22 

2.17 

C7(Io.4o) 

C7(T0.4o) 

7.05 

6.08 

0.20 

3.15 

0.20 

4.15 

Cs(Io.4o) 

C8(T0.4o) 

6.78 

5.05 

6.04 

2.09 

6.04 

6.05 

Cg(Io.4o) 

Cg(To.4o) 

3.05 

1.02 

9.98 

1.10 

5.06 

10.01 

Cio(Io.4o) 

Cio(To.4o) 

12.20 

14.01 

-1.97 

0.99 

-3.93 

3.04 

Cn  (I0.40) 

Ch(T0.4o) 

4.80 

3.19 

1.01 

5.98 

3.91 

0.72 

Table  9:  This  table  illustrates  the  registration  results  upon  registering  the  kth  multiscale  component 
C,(T„,4o)  of  the  noisy  translated  image  T  with  the  kth  multiscale  component  Cfc(Io.4o)  of  the  noisy  image 
I0.40  obtained  via  the  multiscale  decomposition  discussed  in  Section  6.  Here,  we  use  m  =  12  hierarchical 
steps  to  decompose  the  noisy  image,  so  we  perform  m  —  12  registration  simulations.  The  transformation 
parameters  ipx  and  ip y  are  the  X-  and  Y-translation  parameters  of  the  optimal  transformation  ip  produced 
by  the  component-wise  multiscale  registration  algorithm.  Recall  that  the  actual  translation  values  are  13 
in  X  and  17  in  Y. 
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Mutual  information 

Mean  squares 

Normalized  correlation 

4'x 

13.03 

13.03 

13.04 

4/Y 

16.92 

16.57 

16.59 

Table  10:  This  table  gives  the  translation  parameters  'I'x  and  TY  obtained  by  computing  the  weighted 
average  of  the  multiscale  registration  transformations  ipk  for  k  =  0,1,...  11,  as  described  in  Equation 
(7.4).  Recall  that  the  actual  translation  values  are  13  in  X  and  17  in  Y. 


squares,  we  set  the  weights  bk  =  1  for  k  —  0, 1  and  6*.  =  0  for  k  >  2.  Finally,  for  normalized  correlation, 
we  set  the  weights  bk  =  1  for  k  =  0, 1,2,3  and  bk  =  0  for  k  >  4.  In  Table  10,  we  present  the  X-  and 
Y-translation  values  corresponding  to  the  weighted  average  T . 

Note  that  since  the  actual  translation  values  are  13  in  X  and  17  in  Y,  our  multiscale  registration 
techniques  provides  accurate  results  in  the  case  in  which  both  the  fixed  and  moving  images  contain 
significant  levels  of  noise.  Since  both  ordinary  registration  methods  as  well  as  denoising  techniques  failed 
to  produce  acceptable  registration  results,  the  success  of  our  multiscale  technique  is  an  indication  of  its 
general  applicability  and  accuracy,  particularly  for  cases  in  which  other  methods  fail  due  to  the  presence 
of  noise. 
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1.  Introduction 

Radiotherapy  is  an  image-guided  intervention  and  imaging  is  involved  in  every  key  step  of  the 
process,  ranging  from  patient  staging,  simulation,  treatment  planning,  and  radiation  delivery  to 
patient  follow  up  (see  figure  1).  The  evolution  of  radiation  therapy  has  been  strongly  correlated  with 
the  development  of  imaging  techniques.  During  the  early  days  when  Roentgen  first  discovered  x-rays, 
two-dimensional  (2D)  transmission  images  of  the  human  body  provided  unprecedented  imagery  of 
bony  landmarks  which  allowed  radiologists  to  deduce  the  location  of  internal  organs.  Using  planar 
radiographs,  radiologists  planned  cancer  treatments  by  collimating  rectangular  fields  that 
circumscribed  the  presumed  tumor  location.  Additional  blocks  placed  daily  to  match  marks  on  the 
patient’s  skin,  and  later  using  low-temperature  melting  dense  alloys.  The  emergence  of  computed 
tomography  (CT)  in  the  1970s  revolutionized  radiation  therapy  and  allowed  us  to  use  image  data  to 
build  a  3D  patient  model  and  design  3D  conformal  radiation  treatment.  In  general,  3D  conformal 
radiation  therapy  (3D  CRT)  is  a  method  of  irradiating  a  tumor  target  volume  defined  in  a  3D 
anatomical  image  of  the  patient  with  a  set  of  x-ray  beams  individually  shaped  to  conform  to  the  2D 
beam’s  eye  view  (BEV)  projection  of  the  target.  The  reduction  in  normal  tissue  irradiation  when 
moving  from  2D  to  3D  should  theoretically  improve  the  therapeutic  ratio  and  allow  the  tumor  target 
volume  to  be  treated  to  a  higher  dose,  thereby  improving  the  probability  of  tumor  control.  Recent 
technical  advances  in  planning  and  delivering  intensity  modulated  radiation  therapy  (IMRT)  provide 
an  unprecedented  means  for  producing  exquisitely  shaped  radiation  doses  that  closely  conform  to  the 
tumor  dimensions  while  sparing  sensitive  structures  ,_3.  The  development  of  3D  CRT  and  IMRT 
places  more  stringent  requirements  on  the  accuracy  of  beam  targeting.  In  practice,  large  uncertainties 
exist  in  tumor  volume  delineation  and  in  target  localization  due  to  intra-  and  inter-organ  motions.  The 
utility  of  modem  radiation  technologies,  such  as  3D  CRT  and  IMRT,  cannot  be  fully  exploited 
without  eliminating  or  significantly  reducing  these  uncertainties.  The  need  to  improve  targeting  in 
radiation  treatment  has  recently  spurred  a  flood  of  research  activities  in  image-guided  radiation 
therapy  (IGRT). 

While  all  radiation  therapy  procedures  are  image  guided  per  se,  traditionally,  imaging 
technology  has  primarily  been  used  in  producing  3D  scans  of  the  patient’s  anatomy  to  identify  the 
location  of  the  tumor  prior  to  treatment.  The  verification  of  a  treatment  plan  is  typically  done  at  the 
level  of  beam  portals  relative  to  the  patient’s  bony  anatomy  before  patient  treatment.  In  current 
literature,  the  term  of  IGRT  or  IG-IMRT  is  employed  loosely  to  refer  to  newly  emerging  radiation 
planning,  patient  setup  and  delivery  procedures  that  integrate  cutting-edge  image-based  tumor 
definition  methods,  patient  positioning  devices  and/or  radiation  delivery  guiding  tools.  These 
techniques  combine  new  imaging  tools,  which  interface  with  the  radiation  delivery  system  through 
hardware  or  software,  and  state-of-the-art  3D  CRT  or  IMRT,  and  allow  physicians  to  optimize  the 
accuracy  and  precision  of  the  radiotherapy  by  adjusting  the  radiation  beam  based  on  the  true  position 
of  the  target  tumor  and  critical  organs.  With  IGRT,  it  is  also  possible  to  take  tumor  motion  into 
account  during  radiation  therapy  planning  and  treatment.  Because  IGRT  improves  precision,  it  raises 
the  possibility  of  shortening  the  duration  of  radiation  therapy  by  reducing  the  number  of  treatment 
sessions  for  some  forms  of  cancer. 

The  purpose  of  this  article  is  to  highlight  the  recent  developments  of  various  available 
imaging  techniques  and  present  an  overview  of  IGRT.  Stanford  experience  on  various  aspects  of 
clinical  IGRT  will  also  be  presented.  After  reading  this  article,  it  is  hoped  that  the  readers  will  have 
an  overall  picture  of  IGRT  and  find  it  easier  to  navigate  themselves  through  the  subsequent  articles  in 


2 


this  issue,  which  focus  on  providing  technical  details  and/or  specific  clinical  applications  of  the 
available  IGRT  tools. 

2.  Issues  in  IGRT 

In  current  3D  CRT  or  IMRT,  uncertainties  exist  in  many  circumstances,  such  as  tumor  target 
definition,  patient  immobilization  and  patient  breathing  motion,  which  make  it  difficult  to  administer 
a  high  radiation  dose  to  the  planned  location.  The  exact  locations  of  the  boundaries  of  the  tumor 
target  and  the  adjacent  sensitive  structures  are  often  not  known  precisely,  and  a  population-  and 
disease  site-based  safety  margin  is  used  routinely  to  cope  with  a  problem  that  is  otherwise  insoluble. 
An  important  task  of  IGRT  is  to  eliminate  or  significantly  reduce  the  margins  involved  in  defining  the 
clinical  and  planning  target  volume  (CTV  and  PTV,  respectively). 

Many  IGRT  solutions  have  been  proposed  to  attack  various  aspects  of  the  problem.  Briefly, 
IGRT  developments  are  focused  in  four  major  areas:  (1)  biological  imaging  tools  for  better  definition 
of  tumor  volume;  (2)  time-resolved  (4D)  imaging  techniques  for  modeling  the  intra-fraction  organ 
motion;  (3)  on-board  imaging  system  or  imaging  devices  registered  to  the  treatment  machines  for 
inter-fraction  patient  localization;  and  (4)  new  radiation  treatment  planning  and  delivery  schemes 
incorporating  the  information  derived  from  the  new  imaging  techniques.  These  are  discussed  in  more 
detail  in  the  following. 


3.  Tumor  target  volume  definition 

3.1  CT,  MRI,  and  ultrasound  (US)  imaging  techniques 

To  be  able  to  ‘see’  the  extent  of  disease  more  clearly  and  define  the  tumor  target  volume  relative  to 
the  patient’s  anatomy  have  been  among  the  most  important  issues  in  radiation  oncology.  CT  has 
played  a  pivotal  role  in  the  process.  Many  radiation  oncology  departments  have  acquired  dedicated 
CT  scanners.  A  typical  patient’s  3D  CT  data  set  has  more  than  100  axial  slices,  each  of  which 
contains  512x512  pixels.  With  16  bits  per  pixel,  a  CT  data  set  can  easily  run  over  50  megabytes.  CT 
has  many  advantages,  including  high  spatial  integrity,  high  spatial  resolution,  excellent  bony  structure 
depiction,  and  the  ability  to  provide  relative  electron  density  information  used  for  radiation  dose 
calculation.  The  recent  development  of  ultra-fast  multi-slice  CT  has  opened  a  new  dimension  to  CT 
technology  and  allows  time-resolved  (4D)  CT  imaging  of  patient’s  cardiac  and  breathing  cycles. 
Using  array  detectors,  multisection  CT  scanners  can  acquire  multiple  slices  or  sections 
simultaneously  and  thereby  greatly  increase  the  speed  of  CT  image  acquisition.  Currently,  all 
manufactures  are  moving  toward  8-,  16-  and  even  32-slice  CT  technology.  Radiation  oncology 
application  of  4D  CT  will  be  discussed  in  Sec.  4.1. 

MRI  provides  superior  soft  tissue  discrimination,  especially  for  CNS  structures  and  within  the 
abdomen  and  pelvis,  and  has  been  widely  used  in  the  diagnosis  and  tumor  delineation.  MRI  is  also 
utilized  for  virtual  simulation  of  radiation  treatment  for  some  specific  disease  sites.  Physically,  MRI 
involves  the  determination  of  the  bulk  magnetization  of  nuclei  within  a  given  voxel  through  use  of 
radio-frequency  (RF)  radiation  and  magnetic  fields.  In  a  clinical  setting,  MRI  is  typically  employed 
together  with  CT  images  with  the  help  of  image  fusion  software  to  delineate  the  extent  of  the 
malignancy.  As  with  other  imaging  techniques,  MR  technology  has  gone  through  a  series  of 
revolutions  in  the  past  three  decades.  MRI  technology  is  moving  toward  higher  field  strengths  to 
further  improve  the  quality  of  MR  images,  as  evidenced  by  the  installations  of  3T  scanners  in  many 
institutions  (9.4  T  MRI  scanners  has  been  installed  in  a  few  institutions).  Fast  cine  MRI  is  also 
becoming  increasingly  available  and  may  offer  physicians  an  alternative  for  imaging  the  temporal 
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process  of  patient  breathing  or  even  heart  beating.  Figure  2  shows  an  example  of  MRI  images 
acquired  at  two  different  phases  for  a  liver  cancer  patient.  In  addition,  the  development  of  some 
specialized  MRI  scans  has  also  attracted  much  attention.  These  include  diffusion  and  perfusion  MRI, 
dynamic  contrast  MRI,  MR  angiography,  MR  spectroscopic  imaging  (MRSI)  and  functional  MRI 
(fMRI).  The  recent  development  of  diffusion  tensor  imaging  (DTI),  for  instance,  enables  diffusion  to 
be  measured  in  multiple  directions  and  the  fractional  anisotropy  in  each  direction  to  be  calculated  for 
each  voxel.  fMRI  measures  signal  changes  in  the  brain  that  are  due  to  changing  neural  activity.  These 
techniques  enable  researchers  to  make  axonal  and  functional  maps  to  examine  the  structural 
connectivity  of  different  regions  in  the  brain  and  may  allow  better  definition  of  brain  tumors  and 
better  sparing  of  sensitive  regions4. 

Ultrasound  (US)  is  another  useful  imaging  modality  for  radiation  therapy.  US  utilizes  high 
frequency  (1~10  MHz)  sound  waves  to  generate  anatomical  images  which  have  high  spatial 
resolution  and  tissue  characterization  discrimination  power  through  image  texture  analysis.  In 
radiation  therapy,  it  has  been  particularly  useful  in  prostate  imaging.  Transrectal  US  permits  an 
examination/localization  of  the  prostate  gland  5-8  and  is  the  imaging  modality  of  choice  in  guiding  the 
prostate  seed  implant  procedure. 

3.2  Biological  imaging 

Regardless  of  the  course  of  therapy,  current  standard  imaging  modalities  such  as  CT  and  MRI  do  not 
always  provide  an  accurate  picture  of  the  tumor  extent,  especially  in  the  zone  of  infiltration  that  may 
be  the  limiting  factor  in  an  attempt  of  a  radical  treatment  approach.  This  has  been  shown  to  be  the 
case  for  gliomas  before  surgical  intervention  9>  ,0.  It  is  also  true  when  attempting  to  determine  the 
volume  of  residual  tumor  for  additional  therapy  owing  to  problems  in  differentiating  post-therapy 
changes  from  residual  tumor.  Indeed,  the  above-mentioned  imaging  modalities  are  anatomic  in 
nature,  i.e.,  they  provide  snapshot  of  a  patient’s  anatomy  without  biological  information  of  various 
organs  or  structures.  Biological  imaging,  defined  as  the  in  vivo  characterization  and  measurement  of 
biological  processes  at  the  cellular  and  molecular  level,  is  an  emerging  multidisciplinary  field  resulting 
from  the  developments  of  molecular  biology  and  diagnostic  imaging  and  shows  significant  promise  to 
revolutionize  cancer  detection,  staging/re-staging,  treatment  decision-making,  and  assessment  of 
therapeutic  response.  MRSI  and  positron  emission  tomography  (PET)  are  two  valuable  modalities  for 
radiation  therapy  planning.  ’H  MRSI  combines  the  advantages  of  obtaining  biochemical  data  by 
water-suppressed  H  MR  spectroscopy  with  the  spatial  localization  of  that  data.  MR  spectroscopy  is 
useful  in  characterization  of  brain  and  prostate  tumors.  In  the  brain,  for  example,  malignant  tumors 
have  an  increased  rate  of  membrane  turnover  (increased  level  of  choline)  and  a  decreased 
concentration  of  neurons.  Furthermore,  spectroscopy  allows  for  the  non-invasive  monitoring  of  the 
response  of  residual  tumor  to  therapy  and  for  differentiating  tumor  recurrence  from  tissue  necrosis. 
Recently,  Pirzkall  et  al  ’  have  applied  multi-voxel  MRSI  to  assess  the  impact  of  MRSI  on  the 
target  volumes  used  for  radiation  therapy  treatment  planning  for  high-grade  gliomas.  It  was  found 
that,  although  Tj-weighted  MRI  estimated  the  region  at  risk  of  microscopic  disease  as  being  as  much 
as  50%  greater  than  by  MRSI,  metabolically  active  tumor  tissue  still  extended  outside  the  T2  region  in 
88%  of  patients  by  as  much  as  28mm.  In  addition,  T]-weighted  MRI  suggested  a  lesser  volume  and 
different  location  of  active  disease  compared  to  MRSI.  The  discordance  of  high  grade  glioma  target 
volumes  resulting  from  MRI  was  also  observed  in  other  functional  imaging  modalities  such  as  PET 
andSPECT. 

While  there  is  a  growing  body  of  evidence  now  indicating  that  in  vivo  MRSI  provides  unique 
information  on  metabolism  that  will  ultimately  affect  clinical  diagnosis,  choice  and  monitoring  of 
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therapies,  and  treatment  planning,  in  reality,  MRSI  has  been  mainly  remained  a  research  tool 
confined  to  a  small  number  of  academic  institutions13'18.  PET,  on  the  other  hand,  is  more  widely  used 
and  has  been  harnessed  into  the  planning  process  in  many  clinics.  In  general,  PET  has  lower  image 
resolutions  than  CT  images  and,  with  commonly  used  FDG  tracer,  contains  no  anatomic  information 
about  normal  structures.  Information  derived  from  PET  needs  to  be  fused  with  the  corresponding  CT 
images  for  treatment  planning.  The  fusion  of  PET  and  CT  images  are  simplified  with  the  use  of  the 
hybrid  PET/CT  scanner19  20.  Figure  3  shows  the  data  flow  of  a  typical  PET/CT  scanner.  Flybrid 
PET/CT  systems  have  several  positive  features  that  are  absent  in  stand-alone  PET  and  CT  units. 
PET/CT  is  a  hardware-based  image-fusion  technology  that  virtually  eliminates  the  uncertainty  and 
inconvenience  of  currently  available  software  fusion  of  separate  PET  and  CT  images,  which  are  often 
acquired  with  patients  in  different  positions.  It  should  be  emphasized  that  the  PET/CT  unit  is  not 
simply  a  PET  and  CT  combination:  Not  from  the  perspective  of  system  design,  nor  the  practical 
utility.  Other  than  the  fact  that  one  does  not  have  to  go  through  the  cumbersome  and  time  consuming 
software  fusion  process,  it  has  the  advantages  of  simultaneous  availability  of  the  fused  images, 
convenience  to  the  patient  and  the  physician,  increased  physician  confidence  in  interpreting  the  image 
findings,  and  -30%  of  reduction  in  PET  scanning  time  due  to  the  use  of  CT  data  for  PET  attenuation 
correction. 

3.3  Integration  of  biological  imaging  techniques  and  multimodality  image  fusion 

FDG-PET  provides  a  means  to  study  metabolic  activity  of  tumors  in  vivo.  Initial  studies  incorporating 
FDG-PET  into  treatment  planning  have  been  reported  21-24.  Bradley  et  al  23’  24  have  carried  out  a 
prospective  study  to  determine  the  impact  of  functional  imaging  with  FDG-PET  on  target  volumes 
among  non-small  cell  lung  cancer  (NSCLC)  patients  being  considered  for  definitive  radiation 
therapy.  They  found  that  radiation  targeting  with  fused  FDG-PET  and  CT  images  resulted  in 
alterations  in  radiation  therapy  planning  in  over  50%  of  patients  by  comparison  with  CT  targeting. 
The  changes  included  the  alterations  in  the  AJCC  TNM  stage  (31%  of  the  patients  studied)  and 
modification  of  target  volume  (58%  of  the  patients  studied).  In  a  separate  study,  MacManus  et  al.  22 
reported  that  30%  of  patients  with  locally  advanced  NSCLC  became  ineligible  for  curative 
radiotherapy  because  of  detection  of  either  distant  metastatic  disease  or  intrathoracic  disease  too 
extensive  for  radical  radiation.  Recently,  Howard  et  al  25  have  studied  the  value  of  FDG-PET/CT  for 
esophagus  cancer  and  similar  findings  were  reported. 

3.4  Emerging  PET  tracers  for  oncologic  imaging 

While  FDG-PET  has  been  shown  to  be  effective  for  a  number  of  malignancies,  imaging  of  many 
other  neoplasms,  such  as  breast  cancer  and  prostate  cancer,  with  FDG  has  shown  less  success.  26-29 
Many  pitfalls  have  previously  been  described  with  FDG-PET  imaging.  The  FDG  tracer  can  be  non- 
specifically  taken  up  by  several  benign  conditions  such  as  inflammatory  disease,  pneumonia,  brown 
fat,  muscle,  bowel  uptake,  and  granulomatous  disease.  Also,  slow  growing  indolent  tumors  may 
exhibit  only  a  mild  increase  in  glucose  metabolism  and  therefore  be  missed  by  FDG-PET  30"34.  Thus 
FDG-PET  is  only  minimally  useful  for  the  evaluation  of  indolent  tumors  such  as  organ-confined 
prostate  cancer.  The  recent  development  of  [  FJfluorothymidine  (FLT)  ’  provided  a  new 
opportunity  to  improve  the  sensitivity  and  specificity  of  PET  imaging  of  cancer.  Because  there  is 
upregulation  of  thymidine  transport  into  malignant  cells  due  to  accelerated  deoxyribonucleic  acid 
synthesis,  either  nC  or  18F-labeled  thymidine  radiotracers  can  be  used  to  determine  cellular 
proliferation.  Several  studies  have  shown  that  the  accumulation  of  FLT  correlates  better  with 
proliferation  in  comparison  with  the  commonly  used  FDG  tracer  ’  .  Recently,  Smyczek-Gargya 
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et  al.  38  have  reported  FLT-PET  imaging  experiments  involving  12  patients  with  14  primary  breast 
cancer  lesions  (T2-T4).  Thirteen  of  the  14  primary  tumors  demonstrated  focally  increased  FLT 
uptake.  The  study  showed  that  FLT-PET  is  suitable  for  the  diagnosis  of  primary  breast  cancer  and 
locoregional  metastases  and  the  high  image  contrast  of  the  technique  may  facilitate  the  detection  of 
small  foci. 

Agents,  such  as  antisense  molecules,  aptamers,  antibodies,  and  antibody  fragments,  can  be 
aimed  at  molecular  targets  for  biological  imaging.  Tumor  receptors  and  certain  cellular  physiologic 
activities,  including  metabolism,  hypoxia,  proliferation,  apoptosis,  angiogenesis,  and  infection, 
provide  such  targets.  In  addition  to  FLT,  there  are  several  other  new  nuclide  imaging  tracers  under 
clinical  or  laboratory  investigations  30’  39'46,  which  include,  to  name  a  few,  nC-Acetate  47‘50,  l8F- 
choline  51’52,  1 'C-choline53'55,  64Cu-DOTA-[Lys3]Bombesin56,  l8F-FMISO  57'59,  18F-FAZA  60,  64Cu- 
ATSM61.  For  example,  carcinogenesis  is  often  characterized  by  enhanced  cell  proliferation  and 
transformation  and  elevated  levels  of  choline  and  choline  kinase  activity  in  certain  neoplasmic 
dieases  have  motivated  the  development  of  positron-labeled  choline  analogs  for  noninvasive 
detection  of  cancer  using  PET  53, 62 .  Choline  acts  as  a  precursor  for  the  biosynthesis  of  phospholipids, 
e.g.  phosphoatidylcholine,  the  major  components  of  cell  membrane  63.  Several  preliminary  studies 
have  demonstrated  the  potential  of  the  new  tracer  for  prostate  cancer  and  many  other  cancers  49, 53, 55’ 

62,  64, 65 


3.5  Biologically  conformal  radiation  therapy  (BCRT) 

The  current  3D  CRT  or  IMRT  inverse  planning  is  typically  aimed  at  producing  a  homogeneous  target 
dose  under  the  assumption  of  uniform  biology  within  the  target  volume.  In  reality,  it  is  well  known 
that  the  spatial  biology  distribution  (e.g.,  clonogen  density,  radiosensitivity,  tumor  proliferation  rate, 
functional  importance)  in  most  tumors  and  sensitive  structures  is  heterogeneous.  Recent  progress  in 
biological  imaging  is  making  the  mapping  of  this  distribution  increasingly  possible.  This  new 
development  opens  a  new  avenue  of  research,  coined  BCRT  66'70.  The  goal  of  BCRT  is  to  take  the 
inhomogeneous  biological  information  derived  from  biological  imaging  into  account  and  to  produce 
customized  nonuniform  dose  distributions  on  a  patient  specific  basis.  The  simultaneous  integrated 
boost  (SIB)  to  elective  volumes  recently  appearing  in  the  literature  represents  a  simple  example  of 
BCRT. 


To  establish  BCRT,  three  major  aspects  must  be  addressed:  (i)  Determination  of  the 
distribution  of  biological  properties  of  the  tumor  and  critical  structures;  (ii)  Prescription  of  the  desired 
dose  distribution  for  inverse  planning;  and  (iii)  Inverse  planning  to  generate  most  faithfully  the 
prescribed  nonuniform  dose  distribution.  While  the  development  of  molecular  imaging  techniques  is 
critically  important  in  mapping  out  biology  distributions,  the  successful  integration  of  this 
information  into  IMRT  planning  through  steps  (ii)  and  (iii)  is  also  indispensable  to  fully  exploit  the 
obtained  biology  information  to  improve  patient  care.  With  the  optimistic  assumption  that  spatial 
biology  distributions  within  a  patient  can  be  reliably  determined  using  biological  imaging  in  the 
future,  Yang  and  Xing  70  have  established  a  theoretical  framework  to  quantitatively  incorporate  the 
spatial  biology  data  into  IMRT  inverse  planning.  In  order  to  implement  this  method,  they  first  derived 
a  general  formula  for  determining  the  desired  dose  to  each  tumor  voxel  for  a  known  biology 
distribution  of  the  tumor  based  on  a  linear-quadratic  (LQ)  model.  By  maximizing  the  TCP  under  the 
constraint  of  constant  integral  target  dose,  they  obtained 
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where  Z)07  (/)  is  the  desirable  prescription  dose  at  the  voxel  i  with  the  tumor  cell  density, 
radiosensitivity  and  proliferation  rate  given  by  (/>,,  an  yf),  and  Dref  is  the  reference  dose  for  the 
voxel  with  reference  radiobiological  parameters  ( pKf ,  aref,  y rcf ).  For  a  given  disease  site,  the 

radiation  dose  used  in  current  clinical  practice  with  “intent  to  cure”  can  be  used  as  a  good  starting 
point  in  selecting  the  value  of  Dref.  The  relation  is  quite  general  and  can  be  used  as  prescription  dose  to 
guide  an  arbitrary  inverse  planning  objective  function  aimed  at  producing  a  customized  dose  distribution 
in  accordance  with  the  spatial  biology  information. 

4  Intra-fraction  organ  motion:  managing  the  respiratory  motion 

Components  affecting  the  reproducibility  of  target  position  during  and  between  subsequent  fractions 
of  radiation  therapy  include  the  displacement  of  internal  organs  between  fractions  and  internal  organ 
motion  within  a  fraction.  Depending  on  the  disease  site,  these  components  contribute  differently  to 
the  margins  that  are  to  be  added  around  the  CTV  to  ensure  adequate  coverage.  In  the  thorax  and 
abdomen,  intra-fraction  internal  anatomy  motion  due  to  respiration  is  a  principal  cause  for  large 
safety  margins.  Motion  can  distort  target  volumes  and  result  in  positioning  errors  as  different  parts  of 
the  tumor  move  in  and  out  of  the  image  window  with  the  patient’s  breathing  cycle.  Several  studies, 
conducted  to  examine  the  extent  of  diaphragm  excursion  due  to  normal  respiration,  reported  the 
range  of  motion  from  ~0.5  to  4.0  cm  in  the  superioinferior  direction.  As  a  consequence  of  a 
significant  margin  added  around  the  CTV,  a  large  amount  of  normal  tissue  surrounding  the  CTV  is 
irradiated.  Accounting  for  such  motion  during  treatment  has  the  potential  to  reduce  margins,  leading 
to  reduced  radiation  toxicity  and  risk  of  treatment-induced  complications,  and  yielding  room  for  dose 
escalation. 

A  complete  solution  compensating  for  respiratory  motion  should  ideally  start  at  the  simulation 
stage.  There  have  been  several  studies  to  characterize  the  amplitude,  phase  and  periodicity  of  organ 
motion  71-75  using  fluoroscopic  x-rays,  ultrasound  76,  11 ,  and  magnetic  or  RF  markers  78’  79.  The 
development  and  deployment  of  spiral  and  multi-detector  CT  scanners  have  made  practical  the 
acquisition  of  time-resolved  or  4D  CT  images.  The  reconstructed  images  acquired  with  patients  in 
treatment  positions  provide  4D  models  upon  which  geometric  as  well  as  dosimetric  computations  can 
be  performed.  4D  PET  is  also  becoming  clinically  available  80"82.  Treatment-wise,  respiratory  gating 
technology  and  tumor  tracking  techniques  to  synchronize  delivery  of  radiation  with  the  patient's  own 
respiratory  cycle  are  under  intensive  investigations. 

4.1  4D  CT  Imaging 

A  4D  CT  can  be  either  prospective  or  retrospective.  In  the  former  case,  the  scanner  collects  images  at 
only  one  of  the  breathing  phases  of  the  patient  instead  of  scanning  continuously.  The  retrospective  4D 
CT  scan  results  in  multiple  image  sets,  corresponding  to  different  breathing  phases  of  the  patient,  and 
consists  of  three  relatively  orthogonal  processes83'87:  Recording  of  respiratory  signal(s),  acquisition 
of  time-dependent  CT  projection  data,  and  construction  of  a  4D  image  from  these  data.  The  first 
objective  can  be  achieved  by  tracking  a  surrogate  of  respiration-related  organ  and  tumor  motion,  such 
as  tidal  volume  measured  with  a  spirometer  85,  88,  chest  expansion  monitored  by  a  pneumatic  bellows 
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,  or  a  reflecting  external  marker  placed  on  the  abdomen  and  tracked  with  a  camera  .  Time- 
dependent  CT  data  can  be  acquired  by  oversampling  in  either  helical  or  cine  mode,  and  constructing 
several  CT  slices  over  the  full  respiratory  cycle  at  each  axial  location  86, 90.  Finally,  the  respiratory 
signal  and  CT  data  must  be  combined  into  a  4D  series,  providing  a  CT  volume  at  several  points 
throughout  the  respiratory  cycle.  In  this  section,  we  will  focus  primarily  on  the  implementation  of 
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4D  CT  provided  by  the  Varian  Real-time  Position  Management  (RPM)  camera/software  and  the  GE 
Discovery  ST  multislice  PET/CT  scanner. 

4D  CT  patient  setup  proceeds  along  the  same  lines  as  a  standard  3D  CT  exam.  The  patient  is 
immobilized  on  the  scanner  bed,  and  aligned  using  room  and  scanner  lasers.  Sagittal  and  coronal 
scout  images  are  used  to  verify  patient  positioning,  and  the  setup  is  adjusted  as  necessary.  At  this 
stage  of  the  setup,  the  4D  procedure  begins  to  diverge  from  the  3D  exam. 

The  RPM  system  consists  of  an  infrared  source,  CCD  camera,  and  a  reflecting  block.  The  block 
is  attached  to  the  patient’s  abdomen,  typically  just  inferior  to  the  xiphoid  process,  and  the 
anterioposterior  motion  of  the  block  is  captured  by  the  camera.  This  motion  is  analyzed  in  real-time 
by  Varian  software  on  a  computer  connected  to  the  RPM  camera.  The  breathing  pattern  is  recorded 
for  the  duration  of  the  scan,  and  is  referred  to  as  the  "respiratory  trace."  Once  the  scan  has  finished, 
the  software  retrospectively  computes  the  phase  at  each  point  of  the  respiratory  trace  by  determining 
the  location  of  the  peaks  at  end-inspiration,  and  assigning  percentages  to  interpeak  points  based  on  a 
linear  interpolation  of  the  peak-to-peak  distance.  For  example,  under  this  scheme,  end-inspiration 
occurs  at  0%,  while  end-expiration  typically  appears  near  50  -  60%.  The  peak-to-peak  distance  can 
vary  between  respiratory  cycles,  as  can  the  position  of  end-expiration  with  respect  to  end-inspiration. 

Irregularities  in  a  patient’s  respiratory  pattern  can  often  be  reduced  by  encouraging  the  patient 
to  breathe  calmly  and  consistently,  and  then  relying  on  the  patient’s  compliance  during  the  scan.  If 
this  free-breathing  approach  is  insufficient,  the  RPM  software  can  provide  audio  coaching  in  the  form 
of  a  "breathe  in,  breathe  out"  recording,  which  is  manually  or  automatically  timed  to  the  patient's 
natural  rhythm.  Some  groups  have  used  video  feedback  either  alone  or  concurrently  with  audio 
instructions  91 .  While  audio  and  video  coaching  can  help  by  stabilizing  the  respiratory  period, 
amplitude  and  baseline,  they  can  complicate  matters  for  patients  with  compromised  respiratory 
function,  who  find  it  difficult  or  impossible  to  maintain  a  regular  rhythm.  Another  solution  is  active 
breath  control  (ABC)  92"94,  which  uses  modified  ventilator  equipment  to  control  the  airflow,  albeit  at 
the  (possibly  significant)  expense  of  patient  comfort. 

Once  a  sufficiently  regular  breathing  pattern  has  been  established,  the  CT  data  is  acquired  in 
“cine”  mode.  This  is  a  step-and-shoot  technique,  whereby  the  gantry  rotates  several  times  at  each  bed 
position  in  order  to  acquire  data  over  the  full  respiratory  cycle.  The  raw  data  is  partitioned  into  bins 
corresponding  to  a  user-selected  time  interval  (typically  less  than  1/1 0th  the  average  cycle),  and  CT 
slices  are  automatically  reconstructed  from  these  bins.  Because  several  respiratory  points  are 
sampled  at  each  bed  position,  a  4D  CT  scan  can  take  several  times  as  long  as  a  corresponding  3D  CT, 
resulting  in  typically  1500  -  3000  CT  slices  for  a  20  -  40  cm  axial  FOV. 

The  respiratory  and  scan  data  are  combined  at  a  separate  computer,  the  Advantage  Workstation 
(AW)  (GE  Medical  Systems),  which  uses  the  respiratory  trace  to  sort  the  oversampled  CT  slices 
according  to  their  phase.  The  AW  does  perform  the  phase  calculations,  but  rather  relies  on  the  phase 
stamp  computed  by  the  RPM  during  the  creation  of  the  respiratory  trace  file.  Missing  phases  for  any 
couch  position  are  replaced  with  their  nearest  neighbor,  providing  a  sorted  image  without  any  phase 
gaps.  The  user  can  navigate  through  the  data  in  each  axial  direction,  similar  to  standard  viewing 
software,  but  can  also  scroll  through  the  respiratory  phases  from  end-inspiration  to  end-expiration. 
Individual  phases  can  be  subsequently  extracted,  or  combined  into  averaged  or  minimum/maximum 
intensity  projections,  and  exported  to  planning  software  in  the  form  of  standard  DICOM  series. 
These  exported  image  series  form  the  basis  of  4D  treatment  planning. 

4.2  Unresolved  issues  in  4D  CT 

The  AW  sorts  the  data  by  phase  rather  than  amplitude.  If  the  breathing  were  perfectly  regular  from 
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cycle  to  cycle,  then  phase-  and  amplitude-based  sorting  would  give  very  similar  results.  The 
problem  arises  when  there  is  variation  in  amplitude,  period,  or  baseline,  or  when  the  onset  of  end- 
expiration  does  not  occur  at  the  same  point  each  cycle.  When  these  inconsistencies  arise,  the  sorted 
CT  images  may  contain  mismatch  artifacts  at  the  interface  between  bed  positions  (see  Fig.  4).  Recent 
studies  have  investigated  amplitude-based  binning  as  an  alternative  to  the  phase-based  approach,  and 
it  appears  that  amplitude  sorting  can  improve  image  quality  in  many  cases  95-9 1 .  Other  researchers 
have  matched  adjacent  CT  slices  without  using  a  respiratory  trace,  by  maximizing  the  continuity  of 
CT  units  integrated  over  regions  of  interest 90.  Yet  another  promising  approach  involves  interpolating 
the  CT  data  continuously  between  end-cycle  peaks  using  deformable  models  98. 

A  second  issue  arises  in  the  correlation  between  external  fiducial  movement  and  tumor/organ 
motion.  Amplitude  ratios  between  fiducial  and  tumor  displacement  may  vary  from  cycle  to  cycle, 
and  thoracic  and  abdominal  points  may  involve  relative  phase  shifts  71,  ".  These  shifts  may  be 
especially  crucial  for  tumors  near  the  lung,  where  hysteresis  is  prevalent.  Finally,  larger  organs  such 
as  the  liver  can  experience  substantial  deformation  during  inspiration  and  expiration,  which  may  not 
be  adequately  captured  by  rigid-body  interpolation  between  points  in  the  respiratory  cycle  10°. 

Finally,  even  if  the  4D  CT  images  have  been  acquired  without  problem,  there  remains  the  issue 
of  reproducibility  at  treatment 101 .  If  treatment  planning  and  delivery  are  based  on  4D  CT,  there  is  an 
implicit  assumption  that  anatomic  motion  during  treatment  will  match  the  tumor  and  organ  motion 
observed  during  setup.  This  assumption  can  be  checked  to  some  degree  through  frequent  gated  or 
breath-hold  portal  imaging  l02.  On  the  other  hand,  it  is  reasonable  to  assume  the  patient  will  relax 
over  time,  so  that  their  breathing  becomes  shallower  or  changes  tempo.  Indeed,  studies  have 
demonstrated  that  some  patients  exhibit  systematic  respiratory  changes  over  a  several-week  course  of 
radiation  therapy,  even  with  visual  and  audio  coaching  '°3.  These  issues  strike  at  the  heart  of  IGRT, 
and  provide  a  fertile  ground  for  research. 

4D  CT  usually  delivers  more  radiation  dose  than  the  standard  3D  CT,  since  multiple  scans  at 
each  couch  position  are  required  in  order  to  provide  the  temporal  information.  We  have  developed  a 
method  to  perform  4D  CT  scans  at  relatively  low  current,  hence  reducing  the  radiation  exposure  of 
patients  87.  To  deal  with  the  increased  statistical  noise  caused  by  the  low  current,  we  proposed  a  novel 
4D  penalized  weighted  least  square  (4D-PWLS)  smoothing  method,  which  can  incorporate  both 
spatial  and  phase  information.  The  4D  images  at  different  phases  are  registered  to  the  same  phase  via 
a  deformable  model,  whereby  a  regularization  term  combining  temporal  and  spatial  neighbors  is 
designed.  The  proposed  method  was  tested  with  phantom  experiment  (see  figure  5  for  an  example) 
and  patient  study,  and  superior  noise  suppression  and  resolution  preservation  were  observed. 

4.3  4D  PET  and  related  issues 

4D  PET  poses  a  problem  distinct  from  4D  CT,  in  that  signal  is  inherently  limited  by  the  tolerable 
patient  dose.  The  result  is  that  any  PET  scan  requires  a  significant  amount  of  time  per  bed  position 
(usually  a  few  minutes)  in  order  to  acquire  sufficient  data  to  produce  a  good  image.  This  limitation 
makes  it  difficult  to  partition  PET  data  with  the  same  time  resolution  possible  in  4D  CT,  but 
nonetheless  acquisition  methods  are  clinically  available  to  obtain  PET  images  at  end-inspiration  or 
end-expiration.  The  most  common  solution  is  to  gate  the  PET  scan  at  the  desired  respiratory  end¬ 
point,  and  reconstruct  a  single  bin  of  gated  data  104'106. 

Patient  setup  proceeds  in  the  same  manner  as  an  ungated  PET  scan,  and  a  CT  image  is  acquired 
for  attenuation  correction  just  prior  to  the  PET.  At  this  point,  the  RPM  system  monitors  patient 
breathing  by  tracking  the  reflecting  block,  and  the  acquisition  trigger  is  set  by  the  user  to  occur  at 
some  given  point  (say,  end-inspiration)  in  the  cycle.  Each  time  the  RPM  camera  determines  that  the 
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reflecting  block  (and,  by  extension,  the  patient’s  respiration)  reaches  this  point  in  the  respiratory 
cycle,  a  trigger  is  sent  to  the  scanner,  and  data  accumulation  is  initiated.  Gated  PET  differs 
fundamentally  from  the  4D  CT  protocol,  by  elevating  the  RPM  system  to  this  active  role  in  data 
acquisition. 

In  gated  mode,  the  user  is  able  to  select  both  the  width  of  the  acquisition  window  and  the 
number  of  sequential  bins  to  record  each  respiratory  cycle.  The  bin  width  directly  affects  image 
quality,  since  the  signal-to-noise  ratio  within  an  image  asymptotically  approaches  the  square  root  of 
the  signal  level  107.  Multiple  bin  acquisition  allows  the  capture  of  the  full  respiratory  cycle  in  several 
bins,  offering  the  possibility  of  retrospectively  sorting  into  two  or  more  respiratory  phases.  Each  time 
the  RPM  trigger  is  received,  data  is  directed  to  the  initial  bin,  and  then  to  the  remaining  bins 
sequentially  until  the  next  trigger.  This  process  continues  for  the  duration  of  the  scan.  Ideally  the 
scan  duration  would  be  chosen  such  that  the  first  bin  (the  respiratory  point  of  most  interest)  would 
accumulate  as  many  data  points  as  a  comparable  ungated  scan  (i.e.,  divide  the  bin  width  by  the  duty 
cycle).  In  reality,  since  this  would  lengthen  the  typical  PET  scan  by  a  factor  of  4  or  5,  practical 
clinical  considerations  may  require  the  gated  scan  to  be  shortened,  with  corresponding  image 
degredation. 

Once  the  scan  has  finished,  it  is  possible  to  associate  each  bin  (beyond  the  first  bin)  with  a 
corresponding  point  in  the  respiratory  cycle.  Since  the  respiratory  trace  is  recorded  by  the  RPM,  it  is 
a  relatively  simple  matter  to  analyze  the  respiratory  motion  offline  and  make  this  correspondence.  It 
is  also  possible  to  retrospectively  combine  multiple  bins  into  a  single  bin,  merging  all  the  data  to 
create  an  effectively  ungated  scan.  However,  these  methods  are  not  yet  available  from  the  vendor  as 
a  clinical  tool,  and  must  be  performed  by  the  user  in  the  context  of  research  efforts.  Once  the  desired 
bin  has  been  selected,  its  data  can  be  reconstructed  using  the  vendor-supplied  filtered  backprojection 
or  OS-EM  algorithms.  The  image  results  can  subsequently  be  exported  to  treatment  planning  systems 
for  review,  similar  to  ungated  PET  series. 

A  salient  point  in  the  PET  reconstruction  process  is  the  specification  of  the  attenuation 
correction  map.  The  current  clinical  design  uses  the  CT  scan  acquired  just  prior  to  the  PET 
specifically  for  this  purpose.  This  attenuation  correction  CT  can  be  an  acquired  during  either  free 
breathing  or  breath-hold.  Some  research  has  indicated  that  PET  reconstructions  can  be  quite  sensitive 
to  distortions  in  the  attenuation  correction  map  108-1 10,  and  investigations  are  ongoing  into  the  use  of 
4D  CT  or  other  models  to  accurately  account  for  attenuation  80,  .  On  the  Varian/GE  system,  this 

requires  selecting  the  appropriate  images  from  the  4D  CT  on  the  AW,  sending  these  series  back  to  the 
scanner,  generating  the  attenuation  correction  maps  for  each  4D  PET  bin,  and  then  reconstructing 
each  bin  separately.  Once  again,  this  is  a  research  solution,  and  not  yet  available  from  the  vendor  for 
clinical  use. 

4.4  Combining  4D  PET  with  4D  CT  and  enhancement  of  the  performance  of  4D  PET  with  post¬ 
acquisition  data  processing 

Once  the  4D  PET  has  been  acquired  (either  a  single  phase,  or  perhaps  several),  it  is  possible  to  create 
a  4D  PET/CT  80.  This  involves  manually  selecting  the  PET  and  CT  images  with  corresponding 
respiratory  phases  (or  amplitudes),  and  fusing  them  on  viewing/planning  software.  We  have  recently 
developed  a  4D-4D  image  registration  algorithm,  which  allows  us  to  automate  the  process.  If  the  CT 
and  PET  scans  are  acquired  with  the  same  patient  position  on  the  same  exam,  then  the  process  is  a 
particularly  simple  hardware-based  registration.  On  the  Eclipse  treatment  planning  system,  for 
example,  two  images  (not  just  PET/CT,  but  other  modalities  as  well)  can  be  automatically  fused  if 
they  share  the  same  DICOM  coordinates.  If  the  DICOM  coordinates  are  not  identical,  the  registration 
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is  more  difficult,  requiring  manual  or  automated  shifts  and  rotations  to  match  anatomical  landmarks 
or  fiducials.  Fusion  may  be  additionally  complicated  by  organ  deformation  (see  Sec.  6).  At 

the  present  time,  PET/CT  hardware  fusion  for  ungated  scans  is  well  established  and  readily  available 
within  the  clinical  setting19,  20.  4D  PET/CT  registration,  however,  remains  primarily  within  the 
research  domain. 

The  major  issue  in  4D  PET  is  the  lack  of  statistics.  Since  the  collected  photons  are  divided 
into  several  frames,  the  quality  of  each  reconstructed  frame  is  decreased  with  increasing  number  of 
frames.  The  increased  noise  in  each  frame  heavily  degrades  the  quantitative  accuracy  of  the  PET 
imaging.  We  have  recently  developed  two  corrective  methods  to  enhance  the  performance  of  4D 
PET.  The  first  method,  coined  “retrospective”  stacking  (RS)  ’  ’  ,  combines  retrospective 

amplitude-based  binning  of  data  acquired  in  small  time  intervals,  with  rigid  or  deformable  image 
registration  methods.  Unlike  gating  techniques,  RS  uses  data  along  the  entire  respiratory  cycle, 
thereby  minimizing  the  need  for  lengthened  scans  while  providing  a  four-dimensional  view  of  the 
region  of  interest81,82.  In  the  second  approach114,  we  reconstruct  each  frame  with  all  acquired  4D  data 
by  incorporating  an  organ  motion  model  derived  from  4D-CT  images  by  modifying  the  well-known 
maximum-likelihood  expectation-maximization  (ML-EM)  algorithm.  During  the  processes  of 
forward-  and  backward-projection  in  the  ML-EM  iterations,  all  projection  data  acquired  at  different 
phases  are  combined  together  to  update  the  emission  map  with  the  aid  of  the  deformable  model,  the 
statistics  are  therefore  greatly  improved.  Both  phantom  and  patient  studies  have  indicated  promising 
potential  of  the  two  methods. 

4.5  Radiation  treatment  planning  based  on  4D  information 

How  to  maximally  utilize  the  time-resolved  image  information  derived  from  4D  CT  or  PET/CT 
represents  one  of  the  challenges  in  IGRT.  In  reality,  the  information  can  be  integrated  into  radiation 
treatment  planning  and  delivery  at  different  levels.  At  the  lowest  level,  the  4D  images  can  be 
employed  to  determine  the  extent  of  tumor  movement  on  a  patient  specific  basis  and  the  information 
can  then  used  to  design  the  CTV  margin  and  the  radiation  portals  to  accommodate  the  motion.  Figure 
6  shows  an  example  of  lung  patient,  in  which  tumor  boundaries  at  three  distinct  respiratory  phases  are 
plotted.  We  have  referred  to  this  type  of  treatment  as  “3.5-dimensional”  radiation  therapy.  The  4D 
information  can  also  be  used  for  guiding  breath-hold  or  gated  radiation  therapy.  There  is  also  strong 
interest  in  using  the  4D  data  to  establish  a  4D  patient  model  and  then  to  carry  out  a  4D  radiation 
therapy  plan.  These  are  the  subjects  of  the  following  two  sub-sections. 

4.6  Breath-hold  and  respiratory  gating 

Various  methods  have  been  worked  out  to  counteract  respiratory  motion  artifacts  in  radiotherapy 
imaging.  Among  them  are  breath-hold,  respiration  gating,  and  4D  or  tumor-tracking  techniques73, 74, 
92, 94, 115.  Breath-hold  techniques  either  actively  or  passively  suspend  the  patient’s  respiration  and  treat 
the  patient  during  this  interval.  Deep  inspiration  breath-hold,  active  breathing  control  (ABC)  (which 
forces  shallow  breathing  and  thereby  ‘freezes’  the  tumor  motion  for  a  small  part  of  the  treatment  time 
92),  and  self-held  breath-hold  are  suitable  for  different  types  of  therapy  targeting  different  cancers. 
Different  types  of  equipment,  such  as  stereotactic  frames,  fiducial  tracking  systems,  timers, 
respirometer,  RPM,  or  interlocks,  may  be  needed  depending  on  the  method  of  breath-hold. 

Respiration-gating  methods  involve  tracking  the  patient’s  natural  breathing  cycle  and 
periodically  turning  the  beam  on  when  the  patient's  respiration  signal  is  in  a  certain  phase  of  the 
breathing  cycle  (generally  end-inhale  or  end-exhale).  The  patient’s  respiration  is  continuously 
monitored  and  the  beam  switches  off  as  the  tumor  moves  out  of  the  target  range.  Gated  radiation 
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therapy  can  offset  some  of  the  motion  but  requires  specific  patient  participation  and  active 
compliance.  In  gated  treatment  it  is  required  that  the  CT  images  used  for  treatment  planning  faithfully 
represent  the  actual  treatment  situation.  While  gated  CT  acquisition  at  the  treatment  respiratory  phase 
is  possible,  our  gating  protocol  proceeds  by  picking  up  the  CT  data  at  an  appropriate  phase  from  the 
patient’s  4D  CT  acquired  using  the  method  described  above.  The  gating  interval  is  typically  centered 
at  end-expiration  because  of  the  increased  reproducibility  at  this  point,  and  spans  20-30%  of  the 
breathing  period  in  order  to  provide  a  reasonable  duty  cycle.  Treatment  plans  are  optimized  for  this 
phase  range  by  planning  on  an  averaged  composite  of  the  scans  within  the  interval,  and  using 
maximum-  and  minimum-intensity  pixel  views  to  incorporate  intra-gate  margins.  The  averaged, 
maximum-intensity  and  minimum-intensity  composites  for  a  lung  patient  are  displayed  in  Figure  7. 

4.7  Tumor  tracking 

Similar  to  the  establishment  of  a  3D  geometric  modeling  based  on  traditional  CT  data,  the  availability 
of  4D  imaging  information  makes  it  possible  to  build  a  patient  specific  4D  model.  Figure  8  shows  the 
4D  models  for  three  different  patients  98.  In  obtaining  the  models,  a  BSpline  deformable  registration 
technique  (see  Sec.  6)  was  used  to  register  different  phases  of  the  4D  CT.  Ideally,  organ  motion 
represented  by  the  4D  model  can  be  incorporated  into  the  radiation  treatment  plan  optimization  to 
overcome  the  adverse  effect  of  respiratory  motion  on  IMRT  delivery116.  A  few  groups  115'120  have 
explored  the  feasibility  of  MLC-based  tumor  tracking.  However,  the  interplay  between  different 
phases  has  been  ignored  during  the  plan  optimization  in  most  of  these  studies.  Webb  has  presented  a 
technique  to  model  the  dosimetric  effect  of  elastic  tissue  movement  when  modulated  beams  are 

1  i  i 

delivered  .  In  general,  the  quadratic  inverse  planning  objective  function  becomes 

i  k  t  i 

where  ^  is  prescribed  dose  for  Ath  structure,  is  the  importance  factor  and  ^(fj)  is  the  calculated 

dose  in  voxel  i  at  time  t,  and  the  summation  over  t  represents  the  integral  dose  to  zth  voxel.  For  4D 
planning  it  is  necessary  to  know  the  path  of  each  material  coordinate  during  the  treatment,  which 
involves  registering  the  voxels  in  different  respiratory  phases.  This  can  be  achieved  by  using  a 
deformable  registration  algorithm.  The  optimization  of  the  above  objective  function  or  alike  115, 122‘ 
126  can  proceed  in  a  similar  fashion  as  conventional  3D  inverse  planning  to  derive  the  optimal 
trajectories  of  the  movements  of  the  MLC  leaves.  An  aperture-based  optimization  127-129  seems  to  be 
more  adequate  for  dealing  with  the  organ  motion116. 

4D  methods  propose  to  track  the  tumor  with  the  radiation  beam  as  the  tumor  moves  during  the 
respiration  cycle.  These  techniques  require  acquisition  of  some  form  of  respiration  signal  (infrared 
reflective  markers,  spirometry,  strain  gauges,  video  tracking  of  chest  outlines  and  fluoroscopic 
tracking  of  implanted  markers  are  some  of  the  techniques  employed  to  date),  which  is  assumed  to  be 
correlated  with  internal  anatomy  motion.  Fluoroscopy  and  the  cine  model  electronic-portal-imaging 
device  (EPID)  have  been  proposed  as  a  means  for  real-time  guidance1 30'135.  While  tumor  tracking 
seems  to  be  the  ultimate  goal  of  4D  radiation  therapy,  the  real  challenge  is  clarifying  whether  the  4D 
model  is  repeatable  at  the  time  of  fractionated  treatments,  and  determining  how  to  correctly 
synchronize  the  MLC  movements  with  the  patient  breathing.  Real-time  imaging  and/or  adaptive 
approaches  will  likely  play  a  role  in  this  aspect  and  the  issue  will  surely  need  more  research  for  many 
years  to  come. 
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5  Inter-fraction  organ  movement 

5.1  Current  techniques  in  dealing  with  inter-fraction  organ  movement 

Uncertainty  in  patient  setup  has  long  been  known  as  a  limiting  factor  to  conformal  radiation  therapy. 
Currently,  the  accuracy  of  patient  setup  is  verified  by  megavoltage  (MV)  radiograph  acquired  with 
either  radiographic  film  or  EPID  136'138.  The  patient’s  bony  landmarks  are  used  to  guide  patient 
alignment.  Poor  soft  tissue  contrast  and  often  unclear  projection  of  the  bony  anatomy  are  major 
problems  of  the  approach.  To  improve  the  situation,  planar  kV  x-ray  imaging  has  been  implemented 
in  a  variety  forms130135.  While  these  systems  show  significantly  increased  contrast  for  bony  structure 
differentiation,  observing  soft-tissue  detail  remains  problematic  and  correction  of  daily  organ  motion 
is  still  challenging.  Attempts  have  been  made  to  use  CT  imaging  to  facilitate  the  patient  setup 
process.  Along  this  line,  the  offline  adaptive-radiation-therapy  (ART)  strategy  139, 140  and  in-room  CT 
approach  141  have  been  studied.  The  former  method  aims  to  partially  compensate  for  organ  motion  by 
carrying  out  multiple  CT  scans  in  consecutive  days  in  the  first  week  of  treatment.  The  image  data  are 
then  employed  to  construct  a  patient  specific  PTV  model  from  the  composite  of  the  CTVs  with 
inclusion  of  statistical  variations  of  the  observed  motions.  While  beneficial,  the  approach  is  hardly  an 
ideal  solution  for  dealing  with  the  inter-fraction  organ  motion.  It  relies  on  establishing  a  statistical 
ensemble  of  all  possible  setup  scenarios  under  a  strong  assumption  that  a  limited  number  of  off-line 
CT  scans  can  adequately  describe  the  inherently  complex,  often  unpredictable  inter-fraction  organ 
motion.  Even  when  it  is  achievable,  the  ART  margin  is  not  optimized  on  a  daily  basis  and  there  is 
still  room  for  further  improvement.  An  integrated  CT/LINAC  combination,  in  which  the  CT  scanner 
is  located  inside  the  radiation  therapy  treatment  room  and  the  same  patient  couch  is  used  for  CT 
scanning  and  treatment  (after  a  1 80-degree  couch  rotation),  should  allow  more  accurate  correction  of 
interfractional  setup  errors.  Some  major  radiotherapy  vendors  provide  options  to  install  a  CT  scanner 
in  the  treatment  room.  The  overall  precision  of  EXaCT  Targeting™  from  Varian  has  been  evaluated 
by  Court  et  al  141 .  However,  the  approach  assumes  a  fixed  relationship  between  the  LINAC  isocenter 
and  the  CT  images  and  relies  heavily  on  the  mechanical  integrity  of  the  two  otherwise  independent 
systems.  Increased  capital  cost  and  prolonged  imaging  and  treatment  are  other  concerns. 

Other  patient  localization  techniques  available  include  ultrasound-based  methods,  video-based 
surface  tracking,  on-board  cone-beam  CT  or  kV  x-ray  imaging,  CyberKnife  and  Tomotherapy,  etc. 
For  prostate  radiation  therapy,  on-line  ultrasound  imaging  has  gained  substantial  interest142'144  but  in 
practice  has  been  found  susceptible  to  subtle  sources  of  error  and  inter-user  variability.  On-board 
CBCT  holds  promise  to  become  a  robust  integrated  on-line  imaging  technology  that  can  yield 
unambiguous  soft-tissue  detail  at  the  time  of  treatment.  Furthermore,  CT  numbers  correlate  directly 
with  electron  density,  thereby  providing  the  potential  for  reconstruction  of  the  actual  dose  delivered 
on  a  daily  basis,  in  addition  to  simple  anatomic  structure  alignment.  The  details  of  emerging  CBCT 
will  be  presented  in  the  next  section.  The  robotic  CyberKnife™  from  Accuray  Inc.  represents  another 
promising  technology.  The  system  has  a  feedback  mechanism  in  which  motion  of  the  CTV, 
determined  through  the  Accutrak  infrared-x-ray-correlated  imaging  system,  can  be  fed  back  to  the 
robot  to  track  the  CTV  145.  However,  while  this  improves  the  duty  cycle,  there  is  a  finite  time  between 
measuring  tumor  position  and  arranging  the  compensation  for  motion.  Helical  tomotherapy  is  an 
alternative  means  of  delivering  IMRT  using  a  device  that  combines  features  of  a  linear  accelerator 
and  a  helical  CT  scanner146.  The  commercial  version,  the  HI-ART  II™,  can  generate  CT  images  using 
the  same  MV  radiation  beam  that  is  used  for  treatment.  Since  the  unit  uses  the  actual  treatment  beam 
as  the  x-ray  source  for  image  acquisition,  no  surrogate  telemetry  systems  are  required  to  register 
image  space  to  treatment  space.  Objective  measures  of  noise,  uniformity,  contrast,  linearity,  and 
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spatial  resolution,  and  comparison  with  that  of  a  commonly  utilized  CT  simulator  have  recently 
performed  by  Meeks  et  al  14  . 

5.2  CBCT for  patient  localization 

CBCT  based  upon  flat-panel  technology  integrated  with  a  medical  linear  accelerator  has  recently 
become  available  from  linac  vendors  for  therapy  guidance.  The  volumetric  images  may  be  used  to 
verify  and  correct  the  planning  patient  setup  in  the  linac  coordinates  by  comparing  with  the  patient 
position  defined  in  treatment  plan.  Both  kV  and  MV  beams  have  been  utilized  for  the  application. 
The  former  typically  consists  of  a  kV-source  and  flat-panel  combination  mounted  on  the  drum  of  a 
medical  accelerator  l48,  with  the  kV  imaging  axis  orthogonal  to  that  of  MV  therapy  beam.  The  fan- 
beam  and  cone-beam  MV  CT  in  clinical  applications  have  been  reported  by  Meeks  et  al  147  and 
Poulliot  et  al  149,  respectively.  It  appears  that  the  MV  images  contain  sufficient  resolution  of  bone 
and  air  cavities  to  register  them  to  structures  in  the  planning  CT  with  millimeter  precision  148  l49. 

Currently,  CBCT  is  primarily  used  for  guiding  the  patient  setup150,  151.  The  procedure  is  not 
much  different  from  the  current  patient  treatment,  other  than  the  fact  that  the  AP/LAT  portal  images 
are  replaced  by  volumetric  data.  In  figure  9  we  show  3D  CBCT  images  of  a  prostate  case  in  one  of 
the  fractional  treatments  along  with  the  patient’s  planning  CT  image.  It  is  seen  that  soft-tissue 
structures  and  boundaries  are  visible  to  varying  degrees  in  the  CBCT  images.  The  patient  has 
implanted  fiducials,  which  show  up  on  both  CBCT  and  planning  CT.  Our  experience  indicates  that 
the  cone  beam  data  can  clearly  reveal  setup  error,  as  well  as  the  anatomical  deformations  and  other 
physiological  changes.  During  the  patient  setup  process,  the  3D  CBCT  images  are  registered  with  the 
planning  CT  data  through  the  use  of  either  manual  or  automated  3D  image  registration  software  that 
calculates  shifts  in  x-,  y-  and  z-directions  (depending  on  the  manufacturer,  rotations  can  also  be 
included).  The  movements  determined  during  the  registration  represent  the  required  setup  corrections 
that  should  be  applied  to  the  patient.  Both  phantom  and  patient  studies  from  our  group  have  shown 
that  the  volumetric  imaging  is  superior  to  the  conventional  MV  or  kV  AP/LAT  patient  setup 
procedure.  We  note  that,  if  only  translational  shifts  are  permissible,  the  level  of  improvement  is 
generally  within  2mm  as  compared  with  kV  AP/LAT  setup  procedure  (2D/2D  match).  However, 
CBCT  can  readily  detect  rotational  errors  which  otherwise  be  missed  by  the  2D/2D  match.  In  figure 
10  we  show  the  localization  image  for  a  head  phantom  with  kV/kV  2D/2D  match  and  3D/3D  match 
(CBCT/planning  CT).  The  latter  approach  was  found  to  be  sensitive  enough  to  identify  a  rotational 
error  as  small  as  2°. 

In  practice,  much  effort  is  needed  to  improve  the  robustness  and  efficiency  of  the  volumetric 
image  registration  process.  Furthermore,  in  order  to  fully  utilize  the  volumetric  data,  a  new  paradigm 
with  seamlessly  integrated  simulation,  planning,  verification,  and  delivery  procedure  is  urgently 
needed.  Until  this  is  realized  clinically,  the  volumetric  imaging  is  nothing  but  an  expensive  extension 
of  the  already  functional  planar  verification  approach.  The  capital  cost  and  other  related  overheads  do 
not  seem  to  justify  the  marginal  benefit  if  the  volumetric  data  is  simply  used  for  determining  the 
patient  shift  in  the  space.  However,  one  should  not  underestimate  the  potential  of  the  volumetric 
imaging  for  the  future  of  radiation  therapy,  as  it  opens  a  new  avenue  (perhaps  the  only  avenue),  for  us 
to  realize  the  planned  dose  distribution  with  high  confidence  in  clinical  settings. 

A  few  groups  are  working  on  deformable  model  based  segmentation  and  patient  setup 
procedures113,  1  °'152.  When  deformable  registration  is  used,  there  are  a  few  options  to  achieve  the 
registration  depending  on  whether  the  primary  aim  is  to  match  soft-tissue,  or  to  align  3D  bony 
structures.  In  figure  11  we  show  a  patient’s  CBCT  and  planning  CT  registration  results  using 
different  registration  schemes.  The  multiple  choices  result  from  the  fact  that  the  dimensionality  of  the 
patient  data  is  much  greater  than  that  in  the  patient  setup  procedure  and  suggest  that  deformable 
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registration  is  not  the  ultimate  solution  to  volumetric  image-guided  radiation  therapy.  Nonetheless, 
the  technique  improves  the  current  method  151  since  it  partially  takes  into  account  organ  deformation 
by  achieving  the  closest  overlay  match  possible  between  the  planning  and  CBCT  data  sets  according 
to  our  clinical  objective,  and  serves  as  a  useful  interim  solution  before  a  better  integrated  approach 
becomes  available. 

5.3  CBCT-based  dose  verification 

Another  important  application  of  on-board  volumetric  imaging  is  verification  of  dose  delivery.  We 
have  recently  evaluated  the  accuracy  of  CBCT-based  dose  calculation  and  examined  if  current  CBCT 
is  suitable  for  the  daily  dose  verification  of  patient  treatment153,  154.  A  CT-calibration  phantom  was 
first  used  to  calibrate  both  conventional  CT  and  CBCT.  CT  and  CBCT  images  of  the  calibration 
phantom,  an  anthropomorphic  phantom  and  two  patients  (a  lung  and  a  prostate  case)  were  then 
acquired  for  this  study.  Our  results  indicated  that  the  imperfect  quality  of  CBCT  images  has  minimal 
impact  (<3%)  on  the  dosimetric  accuracy  when  the  intra-fractional  organ  motion  is  small.  When 
intra-fractional  organ  motion  is  large  and  motion  artifacts  is  severe  (e.g.,  in  the  case  of  lung  cancer), 
the  dosimetric  discrepancy  due  to  the  poor  image  quality  of  current  CBCT  was  found  to  be  clinically 
significant.  Furthermore,  in  the  latter  case,  we  found  that  it  is  possible  to  use  a  deformable 
registration  algorithm  to  map  the  corresponding  electron  density  information  from  planning  CT  to 
CBCT  and  then  to  proceed  with  conventional  dose  calculation. 

5.4  Respiratory  motion  artifacts  in  CBCT 

Superior  to  the  common  approaches  based  on  two  orthogonal  images,  CBCT  can  provide  high- 
resolution  3D  information  of  the  patient  in  the  treatment  position,  and  thus  has  great  potential  for 
improved  target  localization  and  irradiation  dose  verification.  In  reality,  however,  scatter  and  organ 
motion  are  two  major  factors  limiting  the  quality  of  current  CBCT.  When  CBCT  is  used  in  imaging 
thorax  or  upper  abdomen  of  a  patient,  respiration  induced  artifacts  such  as  blurring,  doubling, 
streaking,  and  distortion  are  observed,  which  heavily  degrade  the  image  quality,  and  affect  the  target 
localization  ability,  as  well  as  the  accuracy  of  dose  verification.  These  artifacts  are  much  more  severe 
than  those  found  in  conventional  CT  exams,  in  which  each  rotation  of  the  scan  can  be  completed 
within  a  second.  On  the  contrary,  in  CBCT  scan,  the  gantry  rotation  speed  is  much  slower,  typically 
40  seconds  to  1  minute  for  a  full  360-degree  scan  in  acquiring  the  projection  data,  which  is  more  than 
10  breathing  cycles  for  most  patients.  In  figure  12  we  show  the  influence  of  the  same  motion  on  a 
regular  “fast”  CT  scanner  and  CBCT  for  a  motion  phantom,  where  it  is  clearly  seen  that  the  motion 
artifacts  are  much  greater  than  that  in  a  fast  scanner. 

In  the  last  decade  considerable  effort  has  been  devoted  to  finding  solutions  to  remove  motion 
artifacts  and  to  obtain  time-resolved  medical  images.  Wang  and  Vannier155  presented  a  patient- 
motion  estimation  and  compensation  technique  for  helical  CT  systems.  Willis  and  Bresler  156  cast  the 
motion  artifact  problem  as  a  time-varying  tomography  problem  and  required  special-purpose 
hardware  to  optimally  sample  the  spatially  and  temporally  band-limited  CT  signal  space.  A 
parametric  model  for  the  respiratory  motion  was  used  in  MRJ,  and  the  motion  artifacts  were 
successfully  reduced  by  modifying  the  reconstruction  algorithm  l57.  Crawford  et  al  158  brought  the 
concept  into  CT  imaging,  and  derived  an  exact  reconstruction  formula  for  motion  compensation  for 
CT  scans.  Generally,  motion  correction  algorithms  that  assume  a  motion  model  work  well  when  the 
motion  conforms  to  the  model,  but  have  limited  success  when  it  does  not.  As  described  in  Sec.  4.1, 
4D  CT  has  been  developed  in  radiation  oncology  application  in  order  to  explicitly  account  for  the 
respiratory  motion.  The  4D  CT  can  be  used  to  derive  a  patient-specific  deformation  field  and  then 
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incorporated  into  the  CBCT  filtered-backprojection  (FBP)  image  reconstruction  process  159.  The 
algorithm  was  tested  with  simulations  at  different  settings  corresponding  to  conventional  CT  and 
CBCT  scan  protocols,  with  translational  motion  and  more  complex  motion,  and  with  and  without 
Gaussian  noise.  In  figure  13  we  show  the  result  for  the  motion  phantom  depicted  in  figure  12.  159 
Since  the  motion  model  is  directly  derived  from  the  patient  images,  it  should  be  more  accurate  than 
other  artificial  modeling,  and  therefore  more  efficient  motion  correction  is  expected.  In  addition  to 
this  approach,  Sonke  et  al  160  developed  a  CBCT  procedure  consisting  of  retrospective  sorting  in 
projection  space,  similar  to  that  used  in  4D  CT  (Sec.  4.1).  The  subsets  of  projection  data  are  then 
reconstructed  into  4D  CBCT  dataset.  To  achieve  a  sufficient  temporal  resolution,  however,  this  will 
require  slowing  down  the  gantry  rotation.  The  assumption  of  periodicity  of  the  respiratory  motion  is 
also  necessary.  Zeng  et  al  161  proposed  a  method  to  estimate  the  parameters  of  a  non-rigid,  free- 
breathing  motion  model  from  a  set  of  projections  of  thorax  that  are  acquired  using  a  slow  rotating 
CBCT  scanner. 


6.  Rigid  and  deformable  image  registration 

Development  of  an  effective  image  registration  technique  has  been  one  of  the  most  important 
research  areas.  Depending  on  the  mathematical  nature  of  the  transformation,  image  registration  is 
divided  into  rigid  and  deformable  registrations.  In  rigid  transformations,  it  is  assumed  that  the 
geometry  of  the  object  is  identical  in  the  two  input  images  and  no  distortion  occurs  in  the  image 
acquisition  process.  A  rigid  transformation  consists  of  six  degrees  of  freedom:  three  displacement 
parameters  and  three  rotational  parameters.  Deformable  registration,  on  the  other  hand,  is  more 
complicated  and  entails  the  modeling  of  voxel  dependent  distortion.  Clinically,  the  need  for  a  robust 
image  registration  algorithm  to  compare/fuse  images  representing  the  same  structures  imaged  under 
different  conditions  or  on  different  modalities  is  ever  increasing  because  of  the  extensive  use  of 
multi-modality  imaging  and  the  emergence  of  new  imaging  techniques  and  methods. 

Computer-based  rigid  image  registration  has  gained  widespread  popularity  in  the  last  decade  and 
is  used  in  routine  clinical  practice.  In  this  approach,  the  matching  of  the  two  input  images  is 
formulated  into  an  optimization  problem  and  the  best  registration  of  the  two  images  is  obtained  by 
iteratively  comparing  various  possible  matches  until  no  better  registration  can  be  found.  The  search 
for  the  optimal  match  of  the  two  input  images  is  usually  gauged  by  a  ranking  function  constructed 
based  on  some  physical  considerations.  Depending  on  the  nature  of  the  input  images,  the  formulation 
of  the  problem  can  be  highly  complicated.  Court  and  Dong  162  used  a  rigid  transformation  for  the 
correction  of  tissue  displacement.  A  deformable  procedure  based  on  the  finite  element  model  (FEM), 
in  which  images  are  described  as  blocks  of  elastic  materials  on  which  forces  apply,  was  proposed  by 
Bharath  et  al  163  and  Brock  et  al.  164  In  this  approach,  the  parameters  that  control  the  behavior  of  the 
elastic  material  and  are  responsible  for  the  conversion  of  forces  into  local  deformations  of  the  elastic 
material  are  Young’s  elastic  modulus  and  Poisson’s  ratio.  Although  powerful,  the  model  has  the 
drawback  that  values  of  the  elasticity  and  density  constant  for  various  tissues  are  not  readily  available 
and  have  to  be  found  by  a  trial  and  error  procedure.  The  method  also  relies  on  using  complicated 
software  to  generate  a  FEM  mesh  and  masks  of  the  involved  structures.  Schreibmann  and  Xing  have 
proposed  a  general  narrow-band  approach  for  deformable  registration  1 13.  Depending  on  the  problem, 
modeling  of  individual  voxel  movement  can  also  be  made  using  either  B-splines  98,  thinplate  splines 
165, 166,  0ptjca]  fiow  aigorithms  167,  or  fluid  flow  algorithms  16  .  Spline  interpolation  is  a  relatively 
simple  approach  and  our  experience  with  the  algorithm  suggested  that  the  free-form  registration  is 
stable  and  accurate  for  dealing  with  IGRT  image  registration  problems  169.  An  improvement  to  this 
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method  can  be  achieved  by  using  a  spline  model  with  the  smoothness  of  the  deformation  field 
assured  by  the  interpolation  between  a  grid  of  fixed  control  points.  A  simple  method  along  this  line  is 
to  deduce  the  spline  coefficients  from  a  set  of  user-defined  control  points,  as  was  done  by  Fei  et  al. 
170  and  Lian  et  al  166  in  warping  and  registration  of  MR  volumes.  Coselmon  et  al  17lused  a  similar 
technique  to  study  the  accuracy  of  mutual-information-based  CT  registration  of  the  lung  at  exhale 
and  inhale  respiratory  states. 

To  facilitate  the  computer  decision-making  process,  image  pre-processing  or  user  interaction 
may  be  required,  especially  when  dealing  with  a  deformable  image  registration.  The  use  of 
homologous  anatomic  landmark  pairs  on  the  two  input  images  or  the  control  points  is  an  example  of 
this.  In  reality,  the  user  must  have  a  detailed  understanding  of  the  patient  anatomy  and  the 
characteristics  of  the  two  modalities  in  order  to  accurately  identify  the  control  points  on  both  images. 
The  point  pairs  are  usually  obtained  interactively  with  the  user  repetitively  exploring  the  input  image 
sets  and  each  time  trying  to  locate  a  point  in  both  of  them.  Due  to  the  3D  nature,  the  process  is  rather 
tedious  and  difficult  to  perform.  Schreibmann  and  Xing  172  have  developed  a  general  method  to 
facilitate  the  selection  of  control  points  for  both  rigid  and  deformable  image  registrations.  Instead  of 
relying  on  the  interactive  selection  of  homologous  control  point  pairs  on  both  model  and  reference 
images,  in  the  proposed  approach  the  user  needs  only  to  identify  some  small  control  volumes  on  the 
model  image  in  a  somewhat  arbitrary  fashion.  This  new  way  of  image  registration  eliminates  the  need 
for  the  manual  placement  of  the  homologous  control  points  and  allows  us  to  register  the  two  images 
accurately  and  efficiently.  The  method  was  applied  to  both  rigid  and  non-rigid  image  registration 
problems  and  our  results  indicated  that  the  registration  is  reliable  and  provides  a  valuable  tool  for 
intra-  or  inter-modality  image  registration.  In  figure  14  we  show  the  registration  result  of  a  rectal 
cancer  patient  who  has  undergone  both  CT  and  FLT-PET  scans.  The  increased  robustness  and 
confidence  in  the  registration  and  the  increased  speed  of  calculation,  especially  in  the  case  of  the 
deformable  registration,  are  important  features  of  the  new  technique.  Compared  to  the  manual  rigid 
registration,  this  method  eliminates  the  nuisance  of  the  control  point  pair  selection  and  removes  a 
potential  source  of  error  in  registration.  Compared  to  the  automated  method,  the  technique  is  more 
intuitive  and  robust,  especially  in  the  presence  of  image  artifacts. 


7.  Clinical  experience  with  IGRT 

Clinically  implemented  IGRT  techniques  at  Stanford  include  4D  CT,  4D  PET,  Varian  OBI  (both 
planar  and  CBCT),  gating,  and  Accuray  CyberKnife.  Several  image-guided  clinical  protocols  are 
under  investigation.  4D  CT/PET  information  are  used  in  about  40%  of  the  thorax  and  upper  abdomen 
cases  for  patient  specific  tumor  margin  definition  in  3.5D  radiation  therapy  or  for  treatment  planning 
of  gated  radiation  therapy.  CBCT  is  mainly  applied  for  patient  setup  in  the  treatment  of  head  and 
neck,  prostate  and  other  pelvic  diseases.  For  these  sites,  the  CBCT  image  quality  is  reasonable  to 
visualize  soft  tissues,  but  the  quality  is  generally  notably  inferior  to  that  of  the  state-of-the-art  multi¬ 
slice  fan  beam  CT  scanner.  Scan  truncation  artifacts  because  the  patient  shadow  does  not  fit  on  the 
detector  and/or  organ  motion  often  cause  Hounsfield  unit  calibration  problems.  While  this  does  not 
seem  to  influence  the  image  registration,  the  use  of  CBCT  for  dose  calculation  should  proceed  with 
caution.  Our  initial  experience  indicates  that,  when  compared  with  traditional  CT-based  calculation, 
the  dosimetric  error  is  typically  less  than  3%  for  prostate  or  head  and  neck  cases  but  could  be 
significantly  greater  in  the  thoracic  region.  Comparison  between  cone  beam  data  and  portal  image 
derived  setup  errors  show  only  slight  differences  (<2  mm).  However,  we  should  note  that  the 
differences  are  derived  purely  based  on  the  use  of  manufacturer-provided  image-fusion  software. 
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which  often  emphasizes  the  high  intensity  voxels  in  bony  structures.  The  next  step  is  to  implement 
soft-tissue  based  setup  corrections  clinically.  In  reality,  volumetric  data  contain  much  more 
information  compared  to  planar  images,  and  CBCT  promises  to  be  more  useful  in  the  future  when  it 
is  better  integrated  with  treatment  planning  and  delivery  systems.  An  ideal  integration  would  be  to  use 
volumetric  image-derived  information  to  “tweak”  or  re-optimize  the  treatment  plan.  This  work  is  still  in 
progress  at  Stanford. 

As  another  example  of  IGRT  treatment,  we  describe  our  phase  I  and  II  pancreatic  tumor  dose 
escalation  protocol.  The  aim  is  to  use  CyberKnife  to  target  pancreatic  tumors  more  precisely  and  to 
limit  the  toxicities  associated  with  treatment.  In  a  phase  I  study,  we  treated  patients  with  a  single 
fraction  of  15,  20,  and  25  Gy  to  unresectable  pancreatic  tumors  using  the  Cyberknife  stereotactic 
radiosurgery  (SRS)  system  (Accuray  Inc,  Sunnyvale,  CA)  m.  To  track  tumor  movement,  we  implant 
fiducial  seeds  percutaneously  into  the  pancreatic  tumor.  Using  the  Accuray  Synchrony  platform,  a 
model  in  which  the  position  of  the  internal  fiducials  is  correlated  with  the  patient’s  respiratory  motion 
is  developed.  The  Cyberknife  is  able  to  make  real  time  corrections  to  compensate  for  tumor 
movement  during  respiration.  Prior  to  treatment,  patients  underwent  4D  planning  CT  scans.  Using 
this  dataset,  we  are  able  to  visualize  how  the  pancreatic  tumor  moves/deforms  through  respiration  and 
compensate  for  these  dynamic  changes  174.  Minimal  acute  gastrointestinal  toxicity  was  observed 
even  at  the  highest  dose.  All  patients  who  received  25  Gy  had  no  further  local  progression  of  their 
tumor  until  death.  In  a  follow  up  phase  II  study,  a  cohort  of  19  patients  were  treated  with  45  Gy 
conventionally  fractionated  radiation  therapy  using  IMRT  to  the  pancreas  and  regional  lymph  nodes 
followed  by  a  25  Gy  Cyberknife  stereotactic  radiosurgery  boost  to  the  primary  tumor  175.  An 
excellent  rate  of  local  control  with  this  therapy  was  confirmed.  Because  of  the  rapid  progression  of 
systemic  disease,  we  did  not  observe  a  significant  improvement  in  overall  survival  as  compared  to 
historic  controls.  However,  most  patients  had  a  clinical  benefit  (decreased  pain,  increased  activity) 
and  decreased  serum  tumor  marker  for  pancreatic  cancer  (CA-19-9)  following  therapy.  To  document 
that  SRS  truly  resulted  in  an  anti-tumor  effect,  we  routinely  obtain  FDG-PET/CT  scans  before  and 
after  treatment.  Figure  15  is  an  example  of  one  such  study.  There  was  intense  metabolic  activity  of 
the  pancreatic  tumor  prior  to  therapy  with  a  near  complete  resolution  of  FDG  uptake  in  this  patient  4 
weeks  following  therapy.  The  technological  challenge  for  IGRT  to  minimize  toxicity  in  this  clinical 
scenario  is  the  precision  delivery  of  high  dose  radiotherapy.  This  cannot  be  accomplished  without 
taking  into  account  the  respiratory  associated  motion  of  pancreatic  tumors.  This  movement  takes 
place  in  multiple  planes  and  can  result  in  tumor  displacement  of  up  to  1-2  cm.  Furthermore,  tumor 
and  organ  deformation  during  respiration  must  also  be  compensated  for  during  therapy. 


8.  Summary 

With  the  development  of  IMRT  during  the  1990s,  radiation  therapy  entered  a  new  era.  This  new 
process  of  treatment  planning  and  delivery  shows  significant  potential  for  improving  the  therapeutic 
ratio  and  offers  a  valuable  tool  for  dose  escalation  and/or  radiation  toxicity  reduction.  The  improved 
dose  conformity  and  steep  dose  gradients  necessitate  enhanced  precision  and  accuracy  in  patient 
localization  and  spawn  the  development  of  IGRT,  in  which  various  metabolic  and  anatomical  imaging 
techniques  are  integrated  into  the  radiation  therapy  process.  The  overall  goal  of  IGRT  is  to  target 
tumors  more  accurately  while  better  sparing  the  normal  tissues.  Much  recent  effort  is  focused  on 
removing  the  uncertainty  in  the  definition  of  the  target  volume  and  in  the  determination  of  the 
position  of  mobile  and  often  deformable  organs.  Biological  imaging  described  in  this  article  will 
allow  us  not  only  to  delineate  the  boundary  of  the  tumor  volume  based  on  the  tumors’  biological 
characteristics  but  also  to  map  out  the  biology  distribution  of  the  cancer  cells,  affording  a  significant 
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opportunity  for  BCRT  treatment  in  the  future.  Developments  of  effective  4D  CT/PET  techniques  will 
provide  effective  means  for  us  to  understand  the  temporal  dependence  of  the  involved  structures  and 
design  the  best  possible  strategy  for  targeting  the  moving  tumor.  Integration  of  various  imaging  tools 
for  off-line  and  on-line  application  is  also  of  paramount  importance,  enabling  us  to  ensure  the 
planned  dose  distributions  can  be  realized  in  the  clinical  setting.  With  the  newly  available  IGRT 
tools,  physicians  will  be  able  to  optimize  radiotherapy  accuracy  and  precision  by  adjusting  the 
radiation  beam  based  on  the  actual  positions  of  the  target  tumor  and  critical  organs  during  radiation 
therapy  planning  and  treatment.  We  should  mentioned  that  IGRT  is  still  in  its  infancy  and  many 
technical  issues  remain  to  be  resolved,  such  as  the  establishment  of  a  robust  deformable  registration 
method,  auto-mapping  of  the  contours  outlined  on  the  planning  CT  to  CBCT  or  to  different  phases  of 
4D  CT,  and  management  of  the  sheer  volume  of  acquired  image  sets  (both  4D  CT/PET  and  CBCT). 
However,  it  is  believed  that  much  of  these  technical  hurdles  will  be  resolved  with  time,  and  that 
IGRT  will  become  the  standard  of  practice  in  the  future  through  the  effort  of  researchers  around  the 
world. 
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CT  images 
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Figure  3  A  schematic  drawing  of  the  data  flow  in  a  hybrid  PET/CT. 
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Figure  4:  Example  of  4D  CT  where  respiratory  cycle  irregularities  have  produced  significant  interbed 
mismatches  near  the  base  of  the  lung. 
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Figure  5.  Motion  phantom  study  for  the  4D-PWLS  method  with  the  thorax  phantom.  The  left  and 
middle  columns  are  the  original  phases  obtained  from  the  GE  Advantage  Workstation,  for  100  mA 
and  10  mA,  respectively;  the  right  column  shows  the  10  mA  phases  after  4D-PWLS  processing.  From 
top  to  bottom  are  phase  0%,  20%,  40%,  60%,  80%,  respectively.  The  red  rectangles  represent  the 
selected  ROI  for  calculation  of  SNRs,  each  of  which  contains  5x5x5  voxels.  PWLS  Smoothed 
10-mA  scan  resulted  more  than  two-fold  increase  in  the  SNR  for  every  phase  of  the  periodically 
moving  phantom.  Similar  results  were  obtained  in  a  patient  4D  CT  study. 


Figure  6  Tumor  contours  for  three  breathing  phases.  The  contours  labeled  as  CT20  and  CT  40  were 
produced  by  applying  the  deformation  field  on  the  tumor  contours  delineated  in  CTO. 
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(C) 

Figure  7.  Composite  scans  of  a  4D  CT  lung  patient,  (a)  Average  pixel;  (b)  Maximum-intensity 
pixel;  (c)  Minimum-intensity  pixel.  The  maximum-intensity  pixel  composite  reveals  the  motion 
extent  of  hyperdense  tissue  (e.g.,  lung  tumor),  while  the  minimum-intensity  pixel  view  provides  the 
motion  range  of  hypodense  regions  (e.g.,  lung  air  volume). 
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Figure  8.  (a):  The  BSpline  grid  superimposed  on  lung  contours  (b):  On  each  node,  deformation  is 
represented  by  arrows,  where  arrow  length  is  proportional  to  the  deformation,  (c)  and  (d):  Same 
analysis  for  two  additional  patients. 
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Figure  9  CBCT  (top)  and  planning  CT  (middle)  for  a  prostate  case.  The  fusion  of  the  two  types  of 
CT  images  is  shown  in  the  bottom  panel. 
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Figure  10.  Setup  localization  image  for  head  phantom  with  kV/kV  2D/2D  match  (top)  and  3D/3D  CBCT 
match  (middle).  The  image  shown  in  the  bottom  panel  illustrate  that  the  CBCT  is  a  sensitive  technique 
capable  of  picking  up  a  2°  rotational  miss-match  between  the  planning  CT  and  CBCT. 
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Figure  1 1  Image  registration  of  CBCT  and  planning  CT  based  on  bony  structure  matching,  soft  tissue 
matching  and  deformable  registration. 
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Figure  12  (a)  Motion  phantom  for  CT  and  CBCT  simulation  study.  The  left  circle  moves 

diagonally  with  an  amplitude  of  1.5  cm  and  a  period  of  3.52  sec.  (b)  Simulated  sinograms  and 
their  corresponding  reconstructed  images  with  standard  FBP  algorithm  when  the  circles  are 
stationary,  (c)  and  (d)  show  the  sinograms  and  their  corresponding  reconstructed  images  for 
1  sec/rotation  acquisition  (conventional  CT  scan  speed)  and  40sec/rotation  acquisition  (on-board 
CBCT  scan  speed),  respectively. 
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Figure  13  (a)  Phantom  and  images  reconstructed  with  motion  correction  for  CT  and  CBCT  settings.  The  three 
images  represent  the  reconstructed  image  of  stationary  phantom  (left),  the  conventional  “fast”  CT  (middle),  and 
the  CBCT  (right),  (b)  and  (c)  Horizontal  profiles  through  the  moving  circle  for  the  images  sown  in  the  middle  and 
right  panels  (blue  curves).  For  comparison,  the  profiles  for  the  stationary  phantom  (left  panel)  and  images 
reconstructed  without  motion  artifacts  removal  mechanism  are  also  plotted  (black  and  red  curves,  respectively). 
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Figure  14  Sagittal,  coronal  and  axial  views  of  the  FLT- 
PET  (1st  row)  and  CT  images  (2nd  row).  The 
checkerboard  of  the  CT  and  FLT-PET  images  after 
registration  using  our  new  method  is  shown  in  the  3rd 
row.  A  3D  view  of  the  registration  is  also  presented, 
where  an  excellent  coincidence  is  observed  between  the 
bony  structures  revealed  in  CT  (white)  and  PET  images 
(orange).  The  right  two  panels  of  the  4-th  row  show  the 
convergence  behaviors  of  our  method  and  the 
conventional  method.  Our  method  leads  to  reproducible 
shifts  in  x-,  y-,  and  z-directions,  and  the  conventional 
approach  leads  to  large  variations  in  the  shifts. 
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Figure  15  FDG-Pet  images  of  a  pancreatic  patient  before  and  after  radiation  therapy. 
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